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Abstract 


This dissertation is focused on electrical conduction behaviour in granular systems with the purpose 
of acquiring a fundamental understanding towards applications of granular materials. The electrical 
properties of granular materials are of great significance in many engineering applications, such as 
powder metallurgy processes, optimization of battery electrodes, piezo-resistive composites for tactile 
sensing, and soil water content measurements. Performance in such systems can be largely influenced 
by complex multi-physics interactions arising from microstructures of granular materials. The bulk of 
this dissertation is built on six published or submitted papers. Each chapter is prefaced by an 
introductory section presenting the motivation for the corresponding paper and its context within the 
greater body of work. 


After project background and related previous work introduced in Chapter 1 and 2, respectively, 
Chapter 3 and 4 deal primarily with the contact properties between rough surfaces. The obtained 
information at the interfacial scale serves as an experimental and numerical basis for modelling inter- 
particle contacts in granular media. Specifically, in Chapter 3, the first and second papers show the 
effects of roughness and fractality on the normal contact stiffness of rough surfaces using 
experimental and numerical approaches, respectively. In Chapter 4 containing the third and fourth 
papers, we extend discussions to the correlation of mechanical and electrical properties with scaling 
analyses, based on experimentally measured normal contact stiffness and electrical contact resistance. 


Chapter 5 with the fifth paper presents the effects of network configuration on macroscopic network 
responses focussing on the dielectric universal scaling behaviour. A unified description has been 
obtained for the emergent scaling properties of network responses for random two-phase systems with 
varying topological configurations. The established network is pivotal for connecting the effective 
behaviour at the surface scale to macro scale responses of granular materials. 


In Chapter 6, the final paper shows a physical picture illustrating experimentally observed alternating- 
current universal scaling in conductive granular systems under different stress states. An effective 
numerical approach incorporating inter-particle interaction has been provided to simulate electrical 
responses of granular materials. The combination of the studies from macro-scale phenomena, 
network topologies, and inter-particle properties is presented leading to new _ physics-based 
constitutive models that contain lower scale information. 


This dissertation presents a new comprehensive understanding of conduction behaviour in granular 
materials by means of a physics-based framework combining features containing both experimental 
and numerical information obtained across various length scales, guiding design and optimisation of 
various granular materials. 
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1 Introduction 


Granular materials are ubiquitous not only in the real nature but also in a variety of engineering 
applications. The electrical conduction of granular materials combines the unique properties of 
individual granules and collective properties arising from interconnections between granules. In this 
chapter, the background of electrical conduction in granular materials and related research are 
introduced. 


1.1 Granular materials 


Granular materials are conglomeration of macroscopically discrete solid particles, which are large 
enough such that they are not subject to thermal fluctuations. The lower size limit for grains in 
granular materials can be as small as nanoscale. On the upper size limit, the physics of granular 
materials may be applied to ice floes where the individual grains are icebergs and to asteroid belts of 
the Solar System with individual grains being asteroids (Duran, 2012). 


Granular materials are commercially important in applications as diverse as pharmaceutical industry, 
agriculture, and energy production (Duran, 2012). According to material scientist Patrick Richard, 
"Granular materials are ubiquitous in nature and are the second-most manipulated material (Patrick et 
al., 2005) in industry (the first one is water)". 


Research into granular materials can be traced back at least to Charles-Augustin de Coulomb, whose 
law of friction was originally stated for interactions between grains (Duran, 2012). The evolution of 
the particles follows Newton's equations, with repulsive forces between particles that are non-zero 
only when there is a contact between particles. Although granular materials are very simple to 
describe they exhibit a tremendous amount of complex behaviour, much of which has not yet been 
satisfactorily explained. They behave differently than solids, liquids, and gases which has led many to 
characterize granular materials as a distinct form of matter. In some sense, granular materials do not 
constitute a single phase of matter but have characteristics reminiscent of solids, liquids, or gases 
depending on the average energy per grain. However, in each of these states granular materials also 
exhibit properties which are unique. Granular materials also exhibit a wide range of pattern forming 
behaviour when excited, e.g., vibrated or allowed to flow, injected heat and current. As such granular 
materials under excitation can be thought of as an example of a complex system (Duran, 2012; 
Herrmann et al., 2013). 


1.2 Electrical conduction in granular materials 


The discovery and characterization of electromagnetic waves by Heinrich Hertz at the end of the 
nineteenth century ignited studies of the interactions of these new waves with various types of 
materials. Among them, the electrical (and heat) transport in conductive granular materials presents 
unique properties and have found numerous applications in industry, such as electronic components 
(Beloborodov et al., 2007; Dyre et al., 2009), combustion reactors (Maniatis et al., 2001), fusion 
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reactors (Gan et al., 2010), batteries (Ott, 2015), piezo-resistive composites for tactile sensing (Stass1 
et al., 2014), metallurgy processes and 3D printing (Arifin et al., 2014) etc. The transport behaviour of 
granular materials is contributed by not only the grains themselves, such as their bulk electrical 
resistivity, elasticity, plasticity, and surface structure, but also the way they are arranged in the 
packing featured by the density, network topology, and state of compaction. 


The granularity brings about new physics over a wide range of length scales from macro down to 
nano and atom scales, further extending the already rich list of remarkable effects exhibited by 
disordered systems (Beloborodov et al., 2007; Mehta, 2007). Recent years have seen remarkable 
progress in the development of micro and nano technologies which allow the design and fabrication of 
granular conductors at smaller length scales (Beloborodov et al., 2007; Mowbray and Skolnick, 2005; 
Murray et al., 2000). By altering the size and shape of granules one can regulate the inter-particle 
interaction, e.g., quantum tunnelling effects (Beloborodov et al., 2007). In particular, by tuning the 
microscopic parameters, one can create granular materials varying from relatively good conductors to 
pronounced insulators, controlled by the material properties of the grains, inter-particle interaction, 
degree of disorder and network topography (Beloborodov et al., 2003). This makes granular 
conductors a new class of artificial materials with tuneable electronic properties, therefore, a perfect 
exemplary system for studying the metal-insulator transition and related phenomena, with 
implications for engineering applications. 


Many efforts have been made to understand electrical conductivity of granular materials composed of 
conductive and insulating constituents. The intensive research in this field is motivated not only by 
their important technological promises but by the appeal of dealing with an experimentally accessible 
model system that is governed by tuneable cooperative effects of disorder, electron correlations, and 
network configuration (Beloborodov et al., 2007; Beloborodov et al., 2003; Mowbray and Skolnick, 
2005). 


1.3 Phase transitions of electrical conduction 


1.3.1 Branly effect 


One of the first major findings of electrical conduction in granular materials is the phase transition 
observed by Edouard Branly in 1890 (Branly, 1892). Metal fillings in a non-conductive container tend 
to present a normal state of high electrical resistance due to the naturally existing of oxide layer. In the 
presence of an electrical potential difference or receiving the transient electromagnetic waves 
generated by an electrical spark at some distance, the initially high electrical resistance drops down by 
several orders of magnitude, becoming practically conductors (Dilhac, 2009; Hirlimann, 2007; Zhai et 
al., 2015). It has been early observed that a permanent and unique percolation path is along with the 
drop in resistance of the granular medium (Dorbolo et al., 2003), leading to the understanding that 
welding of the grains was the result of the coherer effect. Because of this welding of the grains the 
effect is not reversible although the electrical path can be easily broken by gently striking the fillings, 
restoring the original large resistance of the granular medium (Vandembroucq et al., 1997). 
Additionally, mechanical pressure allows better contacts between grains that reduce the original 
resistance, reducing the sensitivity to the external excitation (Zhai et al., 2016a; Zhai et al., 2015). 


This effect was incorporated in practical radio detectors insuring the pioneering developments of 
radio-telegraphy (Hirlimann, 2013; Kryzhanovskii, 1992) including the first wireless radio 
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transmission. During the last decades, interests for the Branly effect resumed (Creyssels et al., 2007; 
Creyssels et al., 2008, 2009; Dorbolo et al., 2007; Falcon and Castaing, 1993; Falcon and Castaing, 
2005; Falcon et al., 2004), focusing on the electrical breakdown of oxide layers on grains (Dorbolo et 
al., 2007; Falcon and Castaing, 2005; Houssa et al., 1998; Vandembroucgq et al., 1997), the modified 
tunnelling effect through the metal-oxide/semiconductor-metal junction (Holm, 1967), the attraction 
of grains by molecular or electrostatic forces (Dorbolo et al., 2003; Dorbolo et al., 2007), and local 
welding of micro-contacts by a Joule heating effect (Béquin and Tournat, 2010; Dorbolo et al., 2002, 
2003; Dorbolo et al., 2007). A global process of percolation within the grain assembly also has been 
investigated (Vandembroucgq et al., 1997). Among the phenomena proposed to explain the coherer 
effect, it is easy to show that some have only a secondary contribution. Particularly, the coherer effect 
was observed with a chain of conductive grains, 1.e., one-dimensional structure, at a sufficiently high 
imposed voltage, in a way similar to the action at a distance of a spark or an electromagnetic wave 
(Falcon and Castaing, 2005), demonstrating that the Branly effect is primarily arising from the inter- 
particle contact rather than percolation mechanism with in the whole granular assembly. 


1.3.2 Frequency-dependent electrical responses 


The other important aspect of electrical conduction in granular materials is the frequency-dependent 
electrical responses, from DC (direct current) to AC (alternating current) regimes. The determination 
of AC conductivity, 0, and permittivity, ¢, of granular materials is a complicated problem that 
depends on many parameters: the statistical distribution of the shape and size of the grains, the applied 
force, and the local properties at the contact scale of two grains, that is, the degree of oxidization and 
surface roughness (Béquin and Tournat, 2010; Bourbatache et al., 2012; Creyssels et al., 2007; Dyre 
et al., 2009; Falcon et al., 2004; Zhai et al., 2016b). 


Conductive granular materials consisting of randomly distributed conducting and insulating phases 
can be considered as disordered dielectric systems, such as ceramic composites (Almond and Bowen, 
2004; Li and Schwartz, 2006), ion/electron-conducting glasses (Almond and Bowen, 2004; Dyre et al., 
2009; Schréder and Dyre, 2008), amorphous semiconductors (Dyre et al., 2009; Elliott, 1987), etc. In 
these systems, various parameters govern the electrical, thermal, chemical, and/or mechanical 
properties of components of a system across multiple scales from molecular up to macroscopic 
length-scale. Experimental and computational research efforts are increasingly conducted in order to 
gain insights into how these properties combine across scales to determine overall system 
performance. Particularly, the AC behaviour of these mentioned systems have been extensively 
studied. A similar conductivity-frequency dependence was proposed by Jonscher (1977b) and is 
known as the “Universal Dielectric Response” (Jonscher, 1977a; Zhai et al., 2017). This emergent 
property does not arise directly from any particular physical or chemical properties of the involved 
components, but rather is a consequence of the way components combine (Almond et al., 2013; 
Almond et al., 2011; Creyssels et al., 2008; Mccullen et al., 2009; Murphy et al., 2006). More 
importantly, the conductivity spectra for ionic and electronic conduction in crystalline and amorphous 
systems (Bakkali et al., 2016; Dyre et al., 2009; Roling, 1998) under a wide range of temperature 
conditions demonstrate a single master curve, suggesting the validity of the time-temperature 
superposition principle. 


These dielectric mixtures have been effectively approximated as a random network of resistors and 
capacitors (Almond and Bowen, 2004; Creyssels et al., 2008; Mccullen et al., 2009) with 
representative conductors exhibiting a constant conductance and dielectric components exhibiting a 
variable complex admittance which is directly proportional to an angular frequency. As well as being 
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important in understanding the electrical properties of granular materials, the approach based on 
resistor-capacitor network provides a useful test bed for identifying different types of emergent 
behaviour, finding the critical parameters associated with the observed network responses, 
determining fundamental drives, and addressing the question regarding the universality of complex 
system behaviour. 


1.4 Research objectives 


This dissertation is focused on conduction behaviour in granular materials with the purpose of 
acquiring a fundamental understanding towards applications of granular energy materials. In granular 
materials, the granules are large enough to possess distinct electronic characteristics compared with 
the overall granule aggregates, but sufficiently small to exhibit network responses in mesoscopic scale. 
A constitutive model unifying disparate systems at various length scales for the electrical conduction 
in granular materials is of considerable importance. This work describes the conduction behaviour of 
granular materials by introducing a physics-based framework combining features containing 
experimental and numerical information obtained across various length scales, thus gaining access to 
a comprehensive model that can be used directly for granular system design and optimisation. 


1.5 Thesis structure 


The project background and related previous work are introduced in Chapter 1 and 2, respectively, 
followed with six published or submitted papers as the bulk of this dissertation. Each paper is 
prefaced by an introductory section presenting the motivation for the corresponding paper. The first 
four papers are concentrated on the rough interfacial properties in Chapter 3 and 4. The obtained 
information at surface scale can serve as the experimental and numerical basis for modelling particle 
contact in granular media. The following fifth paper in Chapter 5 presents the influences of the 
configurations of fabricated networks of granular materials by interactions between granules on 
network electrical responses. The established network is pivotal for connecting the effective 
behaviour at the surface scale to parameters of macro scale of granular materials. In the last paper in 
Chapter 6, the combination of the studies from macro phenomena, network topologies, and inter- 
particle properties are presented leading to new physics-based constitutive models that contain lower 
scale information. Finally, Chapter 7 concludes this work by summarising the findings and 
implications, and provides an outlook on future work in contact mechanics and granular physics. 


2 Related Work 


In this chapter, related work on the electrical conduction of granular materials is reviewed. More 
detailed and topical literature review with more details can be found in following chapters. 


2.1 Effective properties of granular media 


The problem of electrical conduction in heterogeneous materials including granular materials is 
mathematically analogous to the problems of thermal conductivity, permittivity and magnetic 
permeability. The study of these topics dates back to early works of James Clerk Maxwell and Lord 
Rayleigh in the late 19"" century (Carson et al., 2005; Dyre and Schrgder, 2000). Since then, studies of 
heterogeneous dielectric materials have been in a continuous state of development. The heterogeneous 
structure is replaced with a hypothetical material which is homogeneous one-component material with 
equivalent conductivity, diffusivity, coefficient of thermal expansion, mechanical properties etc. 
Numerous models have been proposed in the last century for the estimation of effective properties, 
e.g., elastic properties, thermal properties, electrical and fluid transport properties of composite or 
porous materials. Detailed information can be found in the reviews of (Carson et al., 2005; Clerc et al., 
1990; Landauer, 1978; Pietrak and Wisniewski, 2015; Torquato, 2013; Wang and Pan, 2008). 
Regarding the approaches, it can be distinguished between averaging methods derived analytically for 
simplified microstructures (such as effective medium theory and percolation theory), or Monte-Carlo 
methods, determined based on numerically generated microstructures. 


Widely used for the estimation of effective transport properties are the effective medium 
approximations (EMA), which belong to the class of mean-field theories. They are based on the 
mathematical solution for the disturbance of a homogeneous field due to a spherical inclusion as for 
example derived by Maxwell (Maxwell, 1881). Due to their nature, effective medium approximations 
are unable to accurately predict the properties of heterogeneous materials beyond the percolation 
threshold. Maxwell was the first to give analytical expressions for effective conductivity of 
heterogeneous media in his famous work on electricity and magnetism. He considered the problem of 
dilute dispersion of spherical particles of conductivity k, embedded in a continuous matrix of 
conductivity k,,, where thermal interactions between filler particles were ignored (Bird et al., 2007). 
Maxwell’s expression is as follows: 
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where ¢g is the volume fraction of the filler. Based on this model, Rayleigh (Rayleigh, 1892) 
considered a material in the form of spherical inclusions arranged in a simple cubic array, embedded 
in a continuous matrix. Hamilton and Crosser (1962) included effects of different particle shapes. 
Hasselman and Johnson (1987) emphasized that the interfacial properties between filler and matrix 
tend to have a large influence on the overall conduction behaviour thus should be taken into account. 
Base on Eq. (2.1), different approximations can be derived to estimate the effective conductivity of 
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composite materials. Widely used models are the so-called self-consistent approach and differential 
approximation. For the self-consistent approach, the composite at a certain arbitrary mixing ratio 1s 
considered as a homogenized matrix phase with an effective conductivity (Carson et al., 2005). In the 
second approach, the composite material is assumed to be constructed incrementally by introducing 
infinitesimal changes to an already existing material, leading to differential equations (Bruggeman, 
1935). However, the approximation based on the aforementioned approaches yields unsatisfactory 
results for compositions with widely different conductivities, predicting a spurious percolation 
threshold independent of microstructure (Ott, 2015). In general, effective medium approximations fail 
to predict the properties of multiphase media close to their percolation threshold. Recent work 
(Ambrosetti et al., 2010; Bertei and Nicolella, 2011; Pennetta et al., 2002) further incorporated 
percolation phenomenon, which will be discussed in the following section. 


2.2 Percolation theory 


Percolation theory was first developed by Flory (1941) and Stockmayer (1943) to describe how small 
branching molecules react and form very large macro-molecules (Sahimi, 1994). Percolation 
phenomena arise in a vast field of applications: transport behaviour, mechanical properties, glass 
transition, spread of fire and diseases, etc. (Aharony and Stauffer, 2003). The classical percolation 
theory centres on the occurring sharp transition for site and bond percolation network. This transition 
point is called the percolation threshold which is fundamentally important in estimating the effective 
conductivity in multi-phase medium. When increasing the amount of conductive phase, one 
eventually reaches a point, 1.e., percolation threshold, at which chains of conductive components 
begin to appear. Formation of such conductive channels causes a significant increase in the effective 
conductivity of the medium. This effect is visible as a shift from a flat to a steep slope of the effective 
conductivity plotted versus the volume fraction of conductive phase. 


2.2.1 Conductivity near the percolation threshold 


For the description of the electrical percolation, we consider a mixture of a metallic component 
having a conductivity, dg, and a dielectric component with conductivity, d¢, corresponding to the 
occupied sites/bonds and non-occupied sites/bonds, respectively. Particularly, o¢ demonstrates an 
increasing trend with increasing frequency, approaching infinity. We further consider the scenarios at 
low-frequency region (w — 0) and high-frequency region (w — ©) by incoporating conductor- 
capacitor mixture and superconductor-conductor mixture (Aharony and Stauffer, 2003; Sahimi, 1994) 


In the case of a conductor-capacitor mixture we have oc=0. As the percolation threshold, Pr, is 
approached from above (Pp > Pr),, the effective conductvitiy, o,, for low-frequency scenarios follow: 


0, X Op (Pr — Pr), (2.2) 


where {4 1S a positive constant depending on the dimensionality of the network. While the volume 
fraction of conductive phase, Pp, below the percolation threshold, we have no conductivity, as shown 
in the Fig. 1 





Fig.1. Schematics of the percolation for conductor-capacitor mixture at (a) low-frequency, and (b) 
high-frequency. 


In the other case of a superconductor-conductor mixture where d¢ — © with infinitely high frequency. 
As the percolation threshold is approached from below (P. < Pr), the effective conductivity oj, at high 
frequency can be described: 


On, X Oc (Pr — P.)*, (2.3) 


with the conductivity becoming infinite above the percolation threshold (Aharony and Stauffer, 2003; 
Sahimi, 1994). Here s is a possitive constant depending on the dimensionality of the network. 


2.2.2 Finite—size scaling 


The value of the percolation threshold is mathematically defined for an infinitely large system 
(Sahimi, 1994). Percolation in finite systems deserves discussion because practical applications 
usually involve finite systems (Clerc et al., 1990)According to the finite-size scaling theory, the effect 
of finite size of the system can be shown through the effective values of percolation threshold Pr(L) 
for a finite system of linear size L: 


1 
Pr = Pr(L) ox Lv, (2.4) 


where the value of v depends on the system dimensionality (Sahimi, 1994). The finite size of the 
network cause a shift in the percolation threshold. 


2.2.3 Percolation in binary granular mixtures 


The concept of coordination numbers is central to the practical implementation of percolation theory 
(Chen et al., 2009; Gan and Kamlah, 2010; Gan et al., 2010; Ott et al., 2013; Zhai et al., 2017). 
Broadly speaking, the coordination number represents the number of contacts a particluar particle 
makes with its neighbouring particles. For a mixture of conductive and non-conductive grains 
randomly packed, each grain can be regarded as a lattice site and the contacts between spheres are 
bonds. In percolation theory, it is distinguished between occupied and unoccupied sites. This 
correlated to conducting and non-conducting particles in a binary granular mixture. The percolation 
threshold for a number of common regular lattices and for random-packed mono-sized hard spheres 
are shown in Table 1. 


Table 1. Site-percolation threshold and critical-volume fraction for the common three-dimensional 
lattices and for random-packed hard spheres (Powell, 1979; Scher and Zallen, 1970). 


a Critical Critical 

Lattice en volume saan: Percolation 
number density 
fraction fraction 

2D 
Triangular 6 0.50 0.9069 0.4534 
Square 4 0.59 + 0.2 0.7854 0.46 + 0.02 
Kagome - 0.6527 0.6802 0.4440 
Honeycomb > 0.70 + 0.2 0.6046 0.42 + 0.01 
3D 
diamond - 0.43 + 0.03 0.3401 0.146+0.01 
SC 6 0.3140.02 0.5236 0.162 + 0.01 
bcc 8 0.24 + 0.02 0.6802 0.163 + 0.01 
fcc 12 0,195 +0.01 0.7405 0.144 +0.01 
hep 12 0.195 + 0.01 0.7405 0.144 +0.01 
RPHS 6 0.310 0.59 0.183 + 0.003 


For mono-sized packing, Scher and Zallen (1970) suggested that the volume fraction depends only on 
the coordination number but not on the lattice structure (Powell, 1979). The dependence of the 
percolation threshold of a given mono-sized random packing on reciprocal of coordination is given in 
Fig. 3. 
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Fig.3. The percolation threshold as a function of the reciprocal of the coordination. The number 2 and 
3 after the lattice type indicates coordination to the 2" and 3™ nearest neighbours (Powell, 1979). 


For different-sized binary granular mixtures of randomly packed hard spheres, the determination of 
the critical volume fraction at percolation threshold can be seen in pioneering work done by 
Fitzpatrick et al. (1974), Oger et al. (1986), Powell (1979), Suzuki and Oshima (1985), Bouvard and 
Lange (1991), Kuo and Gupta (1995), and Chen et al. (2009), as is shown in Fig. 4. The percolation 
threshold was related to one certain coordination number for a range of the size ratios (Kuo and Gupta, 
1995). The prediction of the coordination number can be obtained analytically from size ratio and 
volume fraction. An overview of the different approaches is given in Bertei and Nicolella (2011). 
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Fig.4. Comparison of experimental and numerical results from literature (Bouvard and Lange, 1991; 
Fitzpatrick et al., 1974; Martin and Bouvard, 2004; Oger et al., 1986; Powell, 1979; Suzuki and 
Oshima, 1985) of dependence of critical volume fraction on the size ratio of conducting sphere with 
respect to the isolating sphere. 


The average coordination numbers, Zg, of multi-phase randomly packed spheres with size distibutions 
following various rules are found close to 6 (Bernal and Mason, 1960; Suzuki and Oshima, 1985). 
The value is widely used and is in good agreement with the experimentally determined coordination 
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number of 6.24 for monosized hard spheres (Pinson et al., 1998). The particle size ratio shows an 
effect on the coordination number (Bertei and Nicolella, 2011; Pinson et al., 1998). Mechanical 
compression can increase the coordination numbers and decrease perolation threshold (Martin and 
Bouvard, 2004). 


2.2.4 Effective conductivity of binary granular mixtures 


Percolation probabilities represent the likelihood that particles are clustered in ways that form 
connected conduction pathways, which is central to the effective functioning of composite structures. 
For the calculation of the percolation probability of a binary sphere mixture of M phases including 
conductive phase k and other non-conductive phases, an empirical formula has been put forward 
(Chen et al., 2009; Costamagna et al., 1998): 


0.4 


p= (: = (ote) (2.5) 


where P is the percolation probability. Here, coordination number, Z,, (the number of contacts a 
particle of species k has to particles of the same species), is proportional to the surface-area fraction, 
0;, which is of all k particles and the overall mean coordination number Z. The cordination number 
Z,% can be derived by 
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(2.6) 


where r is the radius of the spheres. The surface-area fraction, 6;, is defined in terms of the volume 
fractions @; of the k particles: 


Px/Vr 
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The determination of effective conductivity through consideration of connectivity within the phase k 
can be written in the following from: 


Corp = Ky, Oo (Px — Pr)®, (2:9) 


where P;, is the volume fraction of the phase k, dg is the conductivity of pure material of phase k, wis a 
universal exponent that in three dimensions is 2 (Aharony and Stauffer, 2003), K, can be obtained by 
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Kk, = —————. 
' (1 = Pre 


(2.10) 


So that when P; = 1, Oeg¢ = YOq Where y is a dimensionless parameter indicating the differences 
between conductivity of the pure material and the effective conductivity of granular media with 
identical material (Costamagna et al., 1998). 
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2.3 Discrete element method 


Granular materials are often modelled by either a discrete element method (DEM) or a continuum 
approach, in order to represent the constitutive behaviour of the materials. Three-dimensional 
reconstructions of grain assemblies using DEM provide a great deal of quantitative information about 
the complexity of granular media including the individual inter-particle contact information and the 
topological features during deformation. DEM is becoming widely accepted as an effective method of 
addressing engineering problems in granular and discontinuous materials, especially in granular flows, 
powder technology, and rock mechanics (Cundall and Strack, 1979). However, discrete element 
methods are relatively computationally intensive, which limits the length of the simulation and the 
number of particles involved (Radjai and Dubois, 2011). Due to the vast amount of process 
parameters and environmental influence, using 3D simulations to decouple the contributions of 
different variables of microstructural evolution remains challenging (Shearing et al., 2013). The 
electromechanical coupling has been implemented in DEM approach (Bourbatache et al., 2012; 
Renouf and Fillot, 2008), suggesting that the electrical computation depends on the mechanical 
solution. The inter-particle electrical conductance is computed from the contact forces and the 
resulted in contact area. The overall conductance of the granular assembly can be calculated by 
solving the group equations obtained by applying Kirchhoff-current law on each individual 
conductive particle (Creyssels et al., 2009; Ott et al., 2013; Zhai et al., 2017) 


2.4 Frequency-dependent network responses 


Dielectric relaxation in different materials and systems is one of the most intensely researched topics 
in physics in recent decades (Jonscher, 1999). The subject of relaxation covers all types of stress relief 
in solids including dielectric, mechanical, photoconductive, and chemical and so on. It is possible to 
see a certain commonality of the behaviour for different types of relaxation in the recovery of strain 
on removal of stress and it implies therefore a time and/or frequency dependence (Jonscher, 1999). 


2.4.1 Resistor-capacitor networks 


The granular packing can be simplified to a resistor-capacitor (RC) network (Almond et al., 2013; 
Almond and Bowen, 2004; Almond et al., 2006; Bouamrane and Almond, 2003; Creyssels et al., 2008; 
Murphy et al., 2006). The network responses have been demonstrated to be robust with respect to 
microstructural details and component positions. Therefore, all electrical elements including resistor 
and capacitor are suggested to assign identical values across the network without referring to the 
contact force distribution for the purpose of simplification (Almond et al., 2006; Mccullen et al., 
2009). 


The frequency dependence of network responses can be divided into three regions. At low frequencies, 
the AC conductance of the capacitors is much less than that of resistors. The network conductivity 1s 
thus dominated by random percolation paths of resistors across the network. At intermediate 
frequencies, where the conductance of the capacitors and resistors become comparable, all 
components contribute to the network conductivity. In this emergent frequency range, conductance 
rises aS a power function of frequency. The power law exponent was found to be equal to the 
Capacitive component proportion (Almond and Vainas, 1999; Bouamrane and Almond, 2003; Dyre, 
1993; Jonscher, 1977a). At high frequencies, the conductance of the capacitors would far exceed that 
of the resistors and thus these components can be regarded as short circuits. For cases of high 
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proportion of capacitors, the network conductance is again determined by the random distribution of 
all resistors, effectively bound together by shorting capacitors and opening inductors (Clerc et al., 
1990). 


2.4.2 Mixing rule 


The network responses for low and high frequency regions can be analytically estimated based on 
spectral method (Almond et al., 2013; Jonckheere and Luck, 1998) and the averaging approach 
(Almond et al., 2013; Milton, 1980; Murphy et al., 2006). A useful approach to view the net 
admittance for emergent power-law scaling region is known as Lichtenecker’s rule (Lichtenecker and 
Rother, 1931). The scaling properties for simple equivalence rules (Almond and Vainas, 1999) are 
traditionally given as 


N 
Oy = ». a;0;, (2.11) 


and 
N 
O;° = » a;(o;)~*, (2.12) 


corresponding to parallel and series connections, where o is the admittance and a the volume fraction. 
A suitable description of the whole system is a network of connected RC elements, representing each 
inter-granular contact. The two mentioned simple circuit-equivalence methods of connecting 
components, represent two opposite extremes. It 1s also noted that the two equations are of the general 
form for N-component systems, 


N 
ov = >. a;(o;)”. (2.13) 


This suggests a formula that lies somewhere between the limits of equations shown above, with v of 
+1 and —1 for series and parallel connections, respectively. Considering the logarithmic mixing rule 
(Almond et al., 2006; Murphy et al., 2006), we can further simplify the conductivity of a system with 
two distinct types of components, such as multi-phase materials or an RC circuit, to 


log(o) = a,log(o,) + azlog(o2), (2.14) 
C=, '6,-, (2.15) 


where @, + @2 = 1. For an RC circuit network, the individual admittances are 0, = 1/R and o> = 
iwC, corresponding to resistors and capacitors. The net admittance may then be expressed as 


y=(2)' Gwo)™, (2.16) 


where (iwC)% can be further expressed by (wC)% Jcos (“2") + isin (*)| . The real part of the 


admittance is the conductance, 


G = Re(Y). (2.17) 
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Taking the logarithm of both sides leads to 
log(G) = a,log(w) + B, (2.18) 
es 
2 
between conductance and impedance, with respect to frequency for networks of randomly positioned 


where B = - a,log(R) + azlog(C) + log (cos ( ))]. This suggests a power law relationship 


R and C components. The parameter, a, describing the proportions of capacitive components, 
determines the exponents of the power functions. 


2.5 Contact mechanics 


A wide range of relaxation phenomena are strongly associated with interfacial processes, such as in 
metal—insulator, semiconductor-insulator, electrode—electrolyte and similar systems. Specifically, the 
interfacial roughness is fundamentally important in the conductance in studied granular systems, by 
determining the inter-particle conductance that act as the RC network element. 


A variety of electrical contact resistance models have been developed since rough surface contacts 
were first considered by Archard (1957), and further developed by Ciavarella et al. (2004), Kogut and 
Komvopoulos (2004), Jackson and Streator (2006), and Jackson et al. (2015). In experiments, many 
mechanisms of surface structure evolution have been observed during electrical conduction through 
rough interfaces, including dielectric breakdown of oxide layers, localised current-induced welding, 
chemical disorder arising with random composition and oxidation processes in corrosive 
environments (Creyssels et al., 2007; Falcon et al., 2004). The complex multi-scale morphologies 
exhibited by rough surface structures give rise to difficulties in the direct quantitative evaluation of 
real interfacial contact area, and therefore the contact conductance. The interfacial properties and 
contact mechanics employed in existing models of granular systems are typically simplified (Béquin 
and Tournat, 2010; Bourbatache et al., 2013; Bourbatache et al., 2012). 


A single electrical contact between two grains can be represented by a resistor, capacitor or inductor, 
depending on the inter-particle distance and force. Separating conductive spheres can be simply 
considered as pairs of capacitors. For widely used metallic materials (Bourbatache et al., 2012; Radjai 
and Dubois, 2011; Renouf and Fillot, 2008), a thin surface oxide layer acting as a resistive film with 
typical thickness ranging from | to 100 nm (Kikuti et al., 2004; Proust et al., 2015) can be usually 
found. Consequently, these dielectric layers cause contacts to be non-conductive, under conditions of 
low loads, effectively acting as capacitors. When two spheres are compressed with sufficient pressure, 
surface asperities penetrate the oxide layer thus forming small metal-to-metal contact patches. This 
leads to the conduction of electrical current with a high-level resistance, governed by the constriction 
resistance (Greenwood, 1966; Holm, 1967; Mikrajuddin et al., 1999). The contacting interfaces 
between grains tend to exhibit complex geometries and structures, over a wide range of length scales, 
governing physical properties and interfacial phenomena and giving rise to constriction resistance. 
The constriction resistance due to the convergence and divergence of current flow at a single contact, 
results in a high resistance state. When the radii of contacting asperities are comparable or smaller 
than the average electron mean free path, electrons can only travel ballistically across the micro- 
contacts, e.g., tunneling (Fisher and Giaever, 1961; Mikrajuddin et al., 1999). Consequently, tunneling 
through the oxide layer contributes to the conduction at this stage (Kogut and Komvopoulos, 2004). 
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Under increasing compression, the numerous independent contacting spots increase in size merging 
into larger contacts. When the effective radius of a contacting asperity becomes larger than the mean 
free path of electrons, Holm contact (Jackson et al., 2012) will be the dominant transport mechanism, 
resulting in a relative lower resistance. With further increasing inter-particle force, large zones of 
metal-to-metal contact appear, likely establishing metal chains along several contacting grains, 
presenting inductive properties. Noticeably, the inductive phase may not appear for metals with thick 
oxide layers, such as RuO,, even with high levels of compression (Grimaldi et al., 2002), and for fine 
powders (Creyssels et al., 2008) with high surface-to-volume ratios and hence a high proportion of 
oxide, the inductive phase may not be observed, under insufficient pressure. 
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3 Contact Mechanics at Rough Interfaces 


In a range of applications employing granular materials, such as geothermal energy (Lund et al., 
2005), soil probing (Deschamps et al., 2000), sustainable agriculture (Lal, 2009), slope stability 
(Gunzburger et al., 2005), and foundations for high-rise buildings (Jinli, 2006), interfacial 
characteristics at the finest length scales strongly impact overall system performance. Normal contact 
stiffness is one of the most important physical quantities related to the generalized force displacement 
behaviour of rough surfaces in contact. The contact stiffness indicating the evolution of true contact 
area depends on the applied pressure and is of notable importance for the study of systems involving 
the physical interactions of multiple bodies including granular matter, electrode contacts, and thermal 
contacts, where the interface-localized structures govern overall system performance. In this chapter, 
we experimentally and numerically investigate the influence of surface topology on normal contact 
stiffness. 
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3.1 Paper 1: Experimental investigation on contact stiffness 


This first paper in this chapter shows an experimental investigation into the effects of roughness and 
fractality on the normal contact stiffness of rough surfaces. We considered different types of rough 
surfaces altered by polishing and by five surface mechanical treatments using different sized particles. 
Surface structures were characterised by conventional first- and second-order roughness parameters 
and fractal dimension calculated based on multiple methods. We employ the nanoindentation to 
evaluate the normal contact stiffness with flat tips. The results of our original experiments address 
existing controversy in the literature and show that contact stiffness at rough interfaces follows a 
power-law function with respect to normal force over a wide range of applied loads. We have further 
connected the experimentally obtained power-law exponent and amplitude to surface structure 
characteristics. These correlations are useful in establishing and validating predictive models for 
contact problems for a wide range of rough surfaces, including the curved rough contacts. In Chapter 
5 and 6, the obtained correlations between surface structure and interfacial stiffness will be integrated 
into granular packings, which can be regarded as multi-contact systems, to understand stress- 
dependent responses of granular materials. 


This paper has been published in Experimental Mechanics in 2016. I was the primary researcher and 
author of this paper, being supervised by Dr. Yixiang Gan and Dr. Dorian Hanaor. 
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Abstract The effects of roughness and fractality on the nor- 
mal contact stiffness of rough surfaces were investigated by 
considering samples of isotropically roughened aluminium. 
Surface features of samples were altered by polishing and by 
five surface mechanical treatments using different sized parti- 
cles. Surface topology was characterised by interferometry- 
based profilometry and electron microscopy. Subsequently, 
the normal contact stiffness was evaluated through flat- 
tipped diamond nanoindentation tests employing the partial 
unloading method to isolate elastic deformation. Three indent- 
er tips of various sizes were utilised in order to gain results 
across a wide range of stress levels. We focus on establishing 
relationships between interfacial stiffness and roughness de- 
scriptors, combined with the effects of the fractal dimension of 
surfaces over various length scales. The experimental results 
show that the observed contact stiffness is a power-law func- 
tion of the normal force with the exponent of this relationship 
closely correlated to surfaces’ values of fractal dimension, 
yielding corresponding correlation coefficients above 90 %. 
A relatively weak correlation coefficient of 60 % was found 
between the exponent and surfaces’ RMS roughness values. 
The RMS roughness mainly contributes to the magnitude of 
the contact stiffness, when surfaces have similar fractal struc- 
tures at a given loading, with a correlation coefficient of 
—95 %. These findings from this work can be served as the 
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experimental basis for modelling contact stiffness on various 
rough surfaces. 


Keywords Contact mechanics - Contact stiffness - Rough 
surfaces - Fractal dimension - Nanoindentation 


Introduction 


In contact mechanics, surface morphology plays an essential 
role in determining how solids interact with one another, with 
significance in many processes including friction, wear, ther- 
mal and electrical conduction [1—3]. The classic Hertzian con- 
tact theory considers elastic solids with smooth profiles of 
curved surfaces connecting the applied force and indentation 
depth, by assuming a distribution of normal pressure in the 
contact area. However, this assumption neglects topological 
features of natural surfaces and the interaction between indi- 
vidual contact regions [2]. The existence of surface roughness 
results in only a certain number of peaks or asperities being in 
contact. Since the true contact area at an individual asperity 
can be of nanoscale dimensions, the real contact area between 
two surfaces is typically orders of magnitude smaller than the 
apparent or nominal contact area. As early as the 1950s, 
Bowden and Tabor had already recognized the significance 
of the surface roughness of contacting bodies in contact me- 
chanics. Since then, the mechanics of rough interfaces have 
been broadly studied on the basis of early contributions by 
Archard [4] and Greenwood and Williamson [5]. In the past 
few years, many studies have been carried out to interpret 
parameters of contact area and normal contact stiffness under 
various loading forces, for different surface characteristics 
[6-9]. It has been theoretically and numerically found that in 
non-adhesive contact, the contact area between rough elastic 
surfaces depends linearly on the normal force and 1s inversely 
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proportional to elastic modulus and the root mean square 
(RMS) slope of the surface [10—12]. 

The true interfacial contact area along with roughness pa- 
rameters could, at some level, be utilised to interpret the me- 
chanics of rough surfaces subjected to an applied normal load 
[13, 14]. However, the determination of the real contact area 
between contacting bodies through experimental measure- 
ment is challenging because most natural surfaces exhibit fea- 
tures across a wide range of length scales [15, 16]. Rough 
surfaces possess complex morphologies and geometries that 
are difficult to quantify definitively. 

Numerical analyses, using molecular dynamics, boundary 
element methods, etc., have been conducted and generally 
support the proportionality between normal force and contact 
stiffness, as was first mentioned in the Greenwood- 
Williamson model [5, 8, 10, 17]. However, other studies have 
reported that, for small to medium loads, the logarithm of 
stiffness exhibits close proportionality to the logarithm of the 
applied normal force, i.e. the contact stiffness k is a power 
function of the normal force Fy: k« Fy, which differs from 
the prediction of Greenwood-Williamson model and the work 
mentioned above. This issue has been the subject of contro- 
versy and discussion in the field of contact mechanics for 
some years [7, 18-21]. 

Several related experiments have been conducted on 
multiscale rough surfaces to gain information pertaining to 
interfacial behaviour. Buzio et al. [22] employed atomic force 
microscopy (AFM) to explore the role of surface morphology 
in contact mechanics at the nanoscale. This approach was 
selected owing to its high sensitivity with respect to vertical 
displacements and normal loads. However, observations of 
morphological effects with nanometric tips (for example, 
sharp AFM tips) are encumbered by technical difficulties, 
e.g., uncontrolled long-range adhesive forces dominate incip- 
ient contact while adsorbed water and contaminants may 
smooth atomic corrugation [23, 24]. The digital image corre- 
lation and ultrasound techniques have been employed to mea- 
sure the contact stiffness of real engineering interfaces, reveal- 
ing that contact stiffness depends upon three main factors: 
contact pressure, surface roughness, and surface hardness 
[25]. Relationships between hierarchical surface structures, 
loading conditions and contact mechanics have been studied 
experimentally and computationally using a variety of other 
macroscopic approaches and a range of materials [2 1, 26—32]. 
However, significant difficulties remain in relating parameters 
of contact stiffness to surface structure descriptors. 

The surface morphology of materials is known to be of 
fundamental importance in governing physical properties 
and interfacial phenomena. Consequently, there exists an ever 
increasing number of experimental techniques for the two and 
three dimensional characterisation of rough surfaces along 
with a broad range of numerical descriptors to quantify these 
structures using parameters such as roughness, skewness, 
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kurtosis, curvature, slope, and fractal dimension. Among these 
surface descriptors, the RMS roughness, RMS slope and frac- 
tal dimension are of considerable value in comparative analy- 
sis and quantitative characterisation of surface structures. The 
RMS slope is commonly chosen as a higher order surface 
descriptor, which can be used to interpret surface phenomena. 
For example, optical and tribological properties at rough 1n- 
terfaces have been varyingly correlated with RMS slope [33, 
34]. Moreover, the fractal dimension of surfaces is considered 
to be an important descriptor useful in representing realistic 
three-dimensional surface features over a wide range of length 
scales [35, 36]. The fractal dimension, a cross-scale surface 
descriptor that incorporates localised and macroscopic motr- 
phological surface behaviour, has been attracting increasing 
attention in the characterisation of surfaces and particles 
across a wide range of fields including stereology, powder 
technology, geology, metallurgy, and computer vision [37, 
38]. The advantage of using surface fractality as a cross- 
scale descriptor stems partly from the tendency of first-order 
descriptors (e.g. RMS roughness) to be dominated by the 
highest scale features, while secondary descriptors (e.g. 
RMS slope) tend to be dominated by the finest scale surface 
characteristics. 

In spite of the significant progress that has been made in 
nanoindentation and surface morphology characterisation, 
which facilitate the mechanical and morphological analysis 
of surfaces at micrometre and nanometre scales, the existing 
experimental results of the structure-dependence of contact 
behaviour of surfaces with random multiscale features are 
limited. In this work, we employed flat tip nanoindentation 
to directly observe normal contact mechanics at rough sur- 
faces and, in conjunction with 3D profilometry, established 
relationships between hierarchical surface structures and ex- 
hibited contact stiffness. 


Method 
Sample Preparation and Characterisation 


In this paper, we considered round disk samples made of alu- 
minium alloy 5005 owing to this material’s suitable chemical 
stability and deformability. Three surface treatment methods 
were applied: (1) polishing, (2) abrasive blasting and (3) sur- 
face mechanical attrition treatment (SMAT). Because the alu- 
minium alloy used 1s relatively soft and ductile, the surfaces of 
these samples can be efficiently modified through standard 
polishing treatment, bead blasting and SMAT procedures. 
Samples with polished surfaces were prepared using sequen- 
tial grinding and polishing steps with final polishing using 
1 um diamond suspension. The other two methods employed 
here accomplish physical modification of the surface details at 
different length scales depending on the size of the particles 
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used in the treatment. Both treatments utilised the impact of 
high-speed particles on specimen surfaces. Specifically, the 
operation of bead blasting is propelling forcibly a stream of 
fine glass beads to impact and modify the surface morphology. 
The average sizes of the two selected groups of glass beads 
used in the blasting treatment of the polished surfaces were 50 
and 300 um, respectively. SMAT alters surface features using 
an excitation mechanism to accelerate smooth steel balls and 
project them on the prepared sample surfaces. The surface 
treatments employed will introduce changes in the surface- 
localized microstructure, and thus may alter the material hard- 
ness. However, the material elasticity will not be strongly 
affected by the surface treatments, e.g. SMAT [39-42]. In 
our work we used three different SMAT processes using 1, 2 
and 3 mm diameter balls, respectively. Figure 1 compares the 
samples achieved through the different methods, using a scan- 
ning electron microscope (SEM, Zeiss ULTRA plus), at the 
identical magnification. It can be seen clearly that samples 
after blasting treatment using glass beads of 50 um present 
more complex and rougher texture. 

To further quantitatively analyse the surface morphology, 
the treated aluminium surfaces were scanned using an optical 
surface profilometer (NanoMap 1OOOWLI). A standard 
1024 x 1024 topographical imaging procedure with a vertical 
scanning range of 20 4m was applied with a green LED light 
source (with a wavelength of ~500 nm) to obtain three- 
dimensional topographical maps of the specimen surfaces. 
Subsequently, surface roughness parameters were determined 
from the digitised surface data across multiple scans of 1 x 
1 mm? for each individual sample. The sample surfaces were 
primarily characterised through two roughness descriptors, 
RMS roughness and RMS slope, which are two widely used 
surface descriptors for the characterisation of rough interfaces. 

For the purpose of quantifying and comparing the prepared 
surfaces across a wide range of length scales, a grid of unit 
dimension L is superimposed in order to mesh the obtained 
digitised surface into a number of triangles. The variance of 
surface areas of different samples was calculated at different 
digital resolutions of the scans, as shown in Fig. 2. For exam- 
ple, when L=X/4 with_X being the total scan length, the sur- 
face is covered by 32 triangles of different areas inclined at 
various angles with respect to the projected plane. The areas of 


Fig. 1 SEM images of 
aluminium samples with different 
surface treatments: (a) polishing 
treatment; (b) SMAT using 2 mm- 
sized steel balls; (c) blasted with 
300 um-sized glass beads; (d) 
blasted with 50 um-sized glass 
beads 
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Polished samples 
SMAT (3 mm) 


SMAT (2 mm) 

SMAT (1 mm) 

Glass beads (0.3 mm) 
Glass beads (0.05 mm) 





Fig. 2 The calculated sample surface areas at different scales. The insets 
show the digitised scans used to calculate the surface area for the samples 
that underwent sand blasting using 50 um-sized glass beads 


all triangles are calculated and summed to obtain a total sur- 
face area A, for a given value of L. These triangular elements 
have an equal projected triangle area, although their real areas 
are different and the total surface area As is larger than the 
projected area of the scanned sample 4,,. The grid size is then 
decreased by a successive factor of 2, and the previously men- 
tioned process continues until L corresponds to the distance 
between two adjacent pixel points, 1.e., the highest resolution 
of the digitised data, which is 1 um for the optical profilometer 
employed in the present work. For all of the six sample sur- 
faces, the obtained surface areas A, at various L values exhibit 
an increasing trend as the grid is made finer with a smaller 
value of L. Of the six types of surfaces, the surface area of the 
sample blasted by glass beads shows more evident variation at 
different scales than that of the SMAT prepared samples or 
polished specimen. A power law relationship is found be- 
tween the calculated surface area and the length resolution 
of digitised scan, in the range from 1 to 1000 um. In other 
words, the sample surface structures exhibit self-affinity over 
a certain range of length scales, which can be characterised by 
the surfaces’ fractal dimension. 

In this paper, the concept of fractal dimension provides a 
useful method for representing rough surfaces in terms of 
cross-scales analysis. A number of existing algorithms have 
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been fully developed to determine surface fractal dimension 
from digitised surface data [38, 43]. Here, we used two 
methods to obtain the fractal dimension: (1) Scaled triangula- 
tion; (2) Cube counting. The process of the triangulation 
method based on cross-scale comparison is illustrated in 
Fig. 2, demonstrating the calculated surface areas at various 
scales. The error bars are achieved over five samples for each 
type of surface. The box-counting approach is also applied 
here as a comparison [44, 45]. For both methods, an aniso- 
tropic scaling is applied to normalise three-dimensional scans, 
making both methods suitable to self-affine surfaces. 

As shown in Table 1, six typical surface types with distinct 
surface details were prepared and characterised prior to nano- 
indentation tests. For each surface type, the mean values of 
roughness descriptors were obtained over ten scans and the 
corresponding standard deviations showed these surface mea- 
sures were representative. The resulting values are reported as 
indicated in Table 1 including the mean RMS roughness, 
RMS slope, and fractal dimension. The two calculation 
methods used here to obtain the values of fractal dimension 
are principally in agreement with each other, with the values 
of fractal dimension inverse to the sizes of particles applied in 
SMAT and glass-bead blasting treatments. Moreover, the 
RMS slope and fractal dimension obtained from both methods 
show similar trends, 1.e., smaller particles used to modify the 
surfaces lead to larger RMS slope and fractal dimension 
values. As for the RMS roughness, results show that the pa- 
rameter has no clear correlation with varying particle sizes. 


Contact Stiffness Measurement Using Nanoindentation 


Surface contact stiffness of aluminium samples with different 
surface morphology was analysed by means of nanoindenta- 
tion tests with three flat indenter tips with different diameter 
sizes: 54.1, 108.7, and 502.6 um (SYNTON-MDP, FLT- 
D050, FLT-D100, and FLT-D500), as shown in Fig. 3. The 
tips have sizes comparable with the size of the particles used to 
modify the surface of the specimens and the roll-off wave- 
length found in the power spectra of the surfaces ranging from 
50 to 250 um [43]. The three flat indenter tips were also 


Table 1 Sample surface characterisation with different treatments 

Surface treatment RMS Roughness 
Reus! um 

1. Polish 0.05712+0.00543 1 

2. SMAT3mm Particle size: 3 mm 4.012+0.3674 

3. SMAT2mm Particle size: 2 mm 2.730+0.2554 


1.152+0.3124 
4.179+0.1943 
2,91020;2 159 


4. SMATI1mm Particle size: 1 mm 
5. GB300 um Blasted by glass beads of 300 um 
6. GB50 um Blasted by glass beads of 50 um 
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examined by profilometery, and were found to exhibit an av- 
erage RMS roughness of less than 0.04 um, and thus in the 
following calculations, the tip surfaces were assumed to be 
ideally smooth. The reason for choosing flat tips is that the 
apparent contact area under the tip does not change with re- 
spect to the indentation depth, which is not the case for spher- 
ical or Berkovich tips. All the tips and specimens under tests 
were properly cleaned using water and compressed air to re- 
move any embedded grains and particles. Cleaning with eth- 
anol and heat treatment (around 120 °C) were also applied to 
remove surface contamination and physisorbed moisture. 
When the flat indenter tip first contacts the sample surface, 
the actual contact area is only a small fraction of the nominal 
contact area. The asperities of the sample surface at contact 
regions are squeezed against the flat tip, deforming elastically 
or plastically. For the purpose of evaluating only the elastic 
responses, 10 partial unloading procedures were employed by 
decreasing the applied load by 10 %, in the course of each 
individual test. After each unloading stage, the loading level 
rose by a factor of two to the next unloading stage, until the 
load reached 500 mN where the last unloading was per- 
formed. All the values of stiffness were determined by aver- 
aging over 10 indentation tests at different positions on each 
surface type. Analogous analytical procedures have been 1m- 
plemented experimentally at aluminium interfaces in the past 
using a method of ultrasonic loading [46] and computationally 
in the framework of finite element analyses [47]. The 
unloading stiffness £ (in units of N/m) was defined as the slope 
of the unloading curve, k=dP/dS, where P designates the load 
and S' is the indentation depth, as indicated in Fig. 3. Subse- 
quently, the reduced elastic modulus £,. can be derived from 
the measured unloading stiffness through the relation: 


Jn k 
£, =———}, 
D. 


where A is the nominal contact area of the indenter tip. Equa- 
tion (1) is a fundamental equation for assessing the elastic 
properties in nanoidentation tests, and has been shown to be 
equally applicable in cases of elastic—plastic contact [48—5 1]. 
The applicability of this formula when plastic deformation 


(1) 


RMS slope Rs Fractal dimension Fractal dimension 
Dy / Triangulation Dy / Cube counting 

0.009 101+0.0009640 2.093+0.06176 2.024+0.04101 
0.08280+0.007314 2.185+0.0351 2.149+0.02543 
0.07633+0.00397 2.228+0.01987 2.156+£0.01308 
0.04864+0.01137 2.337+0.01818 2.233+0.01659 
0.2244+0.01537 2.551+0.02170 2.424+0.02572 
0.2023+0.01022 2.626+0.01736 2.351+40.03633 
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FLT-D500 


Aluminium sample 


FLT-D100 
FLT-D050 


Load, P (AU) 





Indentation depth, S (AU) 


Fig. 3. Typical loading-displacement curve of nanoindentation tests on 
the sample surfaces, in arbitrary units (AU). Partial unloading processes 
were applied to isolate elastic contributions to contact stiffness under 
different loading levels. Three flat indenter tips used in the experiments 
(FLT-D050, FLT-D100, FLT-D500) are also illustrated for comparison 


occurs during the indentation experiments has been the sub- 
ject of numerous studies [52, 53]. 

Considering the fact that the measured reduced elastic 
modulus £’. includes contributions from both the tested spec- 
imen and the indenter tip, the contact stiffness, £, (in the unit 
of N/m”), is estimated from 

] l-v.2— 1-v* 


— 2 
Ty Ee bj ( 








where the subscript 7 indicates the property of the indenter tip 
material, the subscript c the properties of the tested specimen 
and v is the material’s Poisson’s ratio. For the diamond indent- 
er tips used in this research, E; and v; are typically 1140 GPa 
and 0.07, respectively. For the aluminium samples employed 
here, v7. is set as 0.3. 


Results and Discussion 
Experimental Data 


Figure 4 presents the typical measured contact stiffness of 
rough surfaces created by sand blasting treatment with 
50 um diameter glass beads (Sample GB50um) obtained from 
nanoindentation tests with three tips of different sizes. The 
stress was calculated from the ratio of loading force to the 
projected area of the indenter tips. Figure 4 shows that, by 
using different sized tips, the stress range extends over several 
orders of magnitude. With the same maximum force 
(~500 mN) provided by the nanoindenter, the maximum stress 
produced with FLT-D050 was around 100 times larger than 
that of the FLT-D500. The stress provided by all the three 
indenter tips ranged from 0.005 to 214.6 MPa, crossing five 
orders of magnitude with the FLT-D500, FLT-D100, FLT- 
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4 Results with -FLT-D050 
@ Results with -FLT-D100 


¢ Results with -FLT-D500 


FLT-D050 


FLT-D100 


- 2.474 | 


LT-D500 





F/A : MPa 


Fig. 4 Contact stiffness, E-, measured from rough surfaces with sand 
blasting treatment of glass beads of 50 um diameter (GB50um). Raw 
data are from multiple indentation tests at different loading levels, F/A, 
where F'is the normal force and A denotes the projected area of the flat tip 


D050 corresponding to the ranges (shown in Fig. 4) of 0.05 
to 2.474 MPa, 0.103 to 52.95 MPa, and 0.419 to 214.6 MPa, 
respectively. The contact stiffness measured over this range of 
applied stresses varies approximately from 0.01 to 55 GPa. 
It can be seen clearly that the measured contact stiffness 
increases with the loading force, for all tested samples, as 
shown in Fig. 5. Although, as mentioned in the introduction 
to the present work, it is the subject of some controversy, the 
contact stiffness has been predicted by certain reports to ex- 
hibit a power-law relation with the applied normal force, 
which can be described as E.« Fy, where Fy is the normal 
force applied on the surface and a is the exponent of the power 
function [13, 19, 21, 22]. At the same applied stress level, the 
surfaces after SMAT (Sample 2-4) or sand blasting treatment 
(Sample 5—6) show a smaller value of contact stiffness with 
respect to that of the polished surface (Sample 1). The surface 
blasted with glass beads of 50 um diameter (Sample 6) 


01.Polish a, =0.4626 
Y 2. SMAT3mm @, =0,5217 
> 3. SMAT2mm @, = 0.5231 
* 4.SMATimm @, =0.5302 
= 5.GB300upm a, = 0.5694 
e 6.GB50pm a, =0.6048 _ad 


Ae dlog(E, / E) 
dlog(F / EA) 





: ° 10° 10° 10° 
F/(EA) 


Fig. 5 Curve fitting for the normalised stiffness, E/E, and normalised 
applied load, F/(F‘A), for six tested surfaces, with F' being the Young’s 
modulus of the tested material, and A the apparent contact area 


SEM 
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presents the lowest contact stiffness of all the six types of 
surfaces. It is also found that the surfaces with the highest 
values of RMS slope and fractal dimension exhibit compara- 
tively lower contact stiffness in the range of applied stress 
levels. 

As the loading force increases towards the high stress re- 
gimes, the slope of the measured contact stiffness decreases. 
The trend is more pronounced as the normal stress approaches 
the bulk compressive strength of the aluminium alloy (shown 
in Fig. 4). This trend matches the expectation that with defor- 
mation of individual asperities the contact stiffness will con- 
verge to a stable value, close to the elastic stiffness of the bulk 
material. 

For each surface, the measured contact stiffness obtained 
through one indenter tip briefly coincides with that of the other 
two indenter tips of different sizes in the overlapping stress 
ranges. However, compared to the results with smaller tips, 
the measured contact stiffness obtained with larger indenter tip 
(FLT-D500) tends to be larger at similar stress levels. A typical 
example is shown in the red dashed circle in Fig. 4. The 
measured contact stiffness at the maximum stress for FLT- 
D500 is evidently larger than the result achieved with FLT- 
D050 under a similar stress level. The indentation of a single 
asperity or fewer asperities, the likely result of the use of the 
50 um tip, occurs with lower apparent stiffness owing to the 
presence of unconstrained neighbouring surface features. 
While the exact nature of deformation mechanisms differs 
between materials it is likely that this phenomenon occurs to 
varying extents across a broad range of systems. 

An alternative explanation to the increased stiffness ob- 
served with the use of larger indenter sizes is the adhesion 
forces between the tip and the rough surface, which are strong- 
ly affected by the true contact area. A larger tip tends to yield a 
relatively higher level of adhesion, resulting in higher bond 
strength and larger apparent contact stiffness. The adhesive 
force can be observed even in the absence of capillary bridges 
if the surfaces are sufficiently smooth [2, 54]. An estimation of 
bond strength can be made at the instant of surface separation 
at the end of the indentation tests, found typically within the 
range of 10 MPa [55, 56]. In our situation, an adhesive inter- 
action between the indenter tip and the testing specimen could 
not be eliminated even after proper cleaning and drying pro- 
cesses, and it can potentially impact the measurements, par- 
ticularly for cases of polished samples subjected to high ap- 
plied loading forces. 

Another factor that is likely to influence the measured con- 
tact stiffness in nanoindentation experiments is the oxide layer 
on the surfaces of aluminium samples. Aluminium alloys 
ubiquitously exhibit a thin passivated oxide layer arising from 
reaction with atmospheric oxygen. This nanoscale oxide layer 
exhibits locally divergent mechanical properties in a region of 
thickness typically less than 10 nm [57, 58]. In this research, 
the penetration depth of the nanoindentation test ranges from 
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1000 to 5000 nm. The smallest depth is observed in one test 
accomplished with FLT-D500 on a polished sample, where 
1025.6 nm is achieved with 500 mN being the maximum 
force. The influence of the oxide layer is thus expected to be 
of limited significance owing to the considerable difference in 
the thickness of the oxide layer and indentation penetration 
depth. 

To compare the contact stiffness of different surfaces, we 
converted the obtained values for contact elastic modulus E.. 
through equations (1) and (2) to non-dimensional values E/E 
by dividing the value of the Young’s Modulus of aluminium 
alloy 5005, E=69.5 GPa. The non-dimensional stress is de- 
fined as F/(EA), where A is the projected area of the corre- 
sponding tip. Note that the fitting curves are achieved, exclud- 
ing the contributions from the measured stiffness under stress 
levels higher than 100 MPa, where the surface shows an ap- 
parent yield phenomenon. The slope a of the fitting curve can 
be defined as the exponent of the power function between the 
normalised contact modulus and normalised stress. For all the 
six sample surfaces, the value of the exponent a varies from 
0.4626 to 0.6048, changing as the fractal dimension and the 
RMS slope increase. As a comparison, this typical value in 
cases of Hertzian contact of two elastic spheres is 1/3. The 
power-law relationship found here experimentally is in good 
agreement with previous theoretical predictions on a quantita- 
tive basis [13, 18, 19, 21]. 


Correlation Analysis 


From the normalised data shown in Fig. 5 the relation between 
the contact modulus E, and applied loading force F' can be 
described as 


2 =0(2). 3) 


The exponent a and amplitude ( are assumed to be constant 
for a given surface structure. Using this form, we can obtain a 
set of coefficients, a and (, corresponding to the tested sur- 
faces from the measurements in Fig. 5. 

Correlations between the exponent a and the evaluated 
surface parameters are described in Fig. 6, where, for all tested 
surfaces, the exponent a is plotted against the aforementioned 
surface descriptors, 1.e., (1) RMS roughness; (2) RMS slope; 
(3) fractal dimension found from a triangulation method, and 
(4) fractal dimension found using a box counting method. 
Specifically, the exponent a achieved across a wide range of 
stress values illustrates a relatively weak correlation to values 
of RMS surface roughness. While in contrast, this exponent 
correlates more closely to values of the fractal dimension with 
the corresponding correlation coefficients above 90 %. In 
Fig. 6(b), a correlation coefficient of around 90 % has also 
been found with the RMS slope, but the exponent a of the 
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Fig. 6 Correlations between the 
stiffness variation exponent a and 
roughness parameters: (a) Non- 
dimensionalised RMS roughness 
Rrus/R’, with R°=1 um; (b) 
RMS slope; (c) Fractal dimension 
obtained with the scaled 
triangulation method; (d) Fractal 
dimension obtained with the cube 
counting method. Horizontal and 
vertical error bars show the 
corresponding standard 
deviations from the surface 
measurement and curve fitting. 
The sample number shown in the 
figures can be referred to Table | 
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Coefficient of correlation: 0.6414 








SMAT samples yields similar values and does not show strong Figure 7 illustrates the significance of RMS roughness, 
correlation to the RMS slope. This, to a certain degree, is | RMS slope and fractal dimension in governing the coefficient 
supported by numerical investigations of these descriptors re- _@ for the contact stiffness measured with flat tipped diamond 


ported in the literature, which, it should be noted, focused indenters. It is found that the RMS roughness dominates (6 
mainly on the evolution of contact area and surface separation — with the coefficient of correlation being —95 %, whilst the 


rather than stiffness [59]. 


Fig. 7 Correlations between the 
coefficient @ and roughness 
parameters: (a) Non- 
dimensionalised RMS roughness 
Rrvs/R°, with R° being 1 um; (b) 
RMS slope; (c) Fractal dimension 
obtained with the scaled 
triangulation method; (d) Fractal 
dimension obtained with the cube 
counting method. Horizontal and 
vertical error bars show the 
corresponding standard 
deviations from the surface 
measurement and curve fitting. 
The sample number shown in the 
figures can be referred to Table | 





RMS slope and the fractal dimension are not as well correlated 
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Table 2 Correlations between the stiffness parameters and surface roughness parameters 


RMS Roughness R*=Rpys/ 1 um RMS slope Rs 


fe’ a=0.019R*+0.487 
B G=—-1.02R*+15.3 


a=0.503Rs5+0.481 
B=-14.1R5+14.27 


with this coefficient. The RMS roughness is a parameter indi- 
cating primarily the vertical scale of the rough surface. If all 
the tested surface structures, along with indentation depths, 
were to be stretched or compressed by a selected factor yield- 
ing identical roughness ranges and the same value of (3, the 
contact behaviour would be mainly controlled by the fractality 
of the surface as described by the value of its fractal dimen- 
sion. The fitting functions for the stiffness parameters, 1.e., a 
and 3, using surface roughness parameters, including RMS 
roughness, RMS slope and fractal dimension are detailed in 
Table 2. 

Similar formulas have been used as the theoretical approx- 
imation for the contact stiffness from Pohrt and Popov [19], 
Jiang et al. [21] and Komvopoulos and Ye [7]. In all the afore- 
mentioned theories, the relationship between the stiffness and 
pressure was found to follow a power law, with the exponent 
strongly correlated to the value of fractal dimension. Pohrt and 
Popov [19] established a power law based on boundary ele- 
ment simulations for a rigid square indenter over a range of 
normal loads spanning four orders of magnitude, concluding 
that the contact stiffness for low to medium loading conditions 
is most appropriately approximated by a power-law depen- 
dence with the exponent ranging from 0.51 to 0.77, corre- 
sponding to a variation of fractal dimension from 2 to 3 
[19]. It is worth noting that the approximated value of the 
exponent a is 0.2567D, which is rather close to the value of 
0.2168D, found experimentally in the current work. Addi- 
tionally, Jiang et al. [21] have presented a fractal model which 
can be applied to analyse the contact stiffness of machined 
plane joint surfaces, in which the exponent a of the power 
law is described as D/(3—D,). Under a given load, smaller 
values of both surface roughness and fractal dimension can 
increase the true contact area, thereby resulting in a larger 
number of micro-contacts yielding elastic deformation, which, 
in turn, increases the normal and tangential contact stiffness. 


Conclusions 


In this work, the effects of surface roughness and fractality on 
the normal contact stiffness were experimentally demonstrat- 
ed for various rough surfaces. A wide range of applied stress 
across five orders of magnitude is achieved using three flat 
indenter tips of various sizes. The results in our experiments 
show that the contact stiffness follows a power-law function 
with respect to the normal force, with the exponent of this 








Fractal dimension Dy, Triangulation 


a=0.217Dp, +0.029 
B=—-4.28Dp) +228 


Fractal dimension Dy Cube counting 


a=0.301D,+0.133 
B=-7.61D 9429.7 


relationship ranging from 0.4626 to 0.6048, corresponding 
to the polished surface and the surface treated with the finest 
beads. This relationship is then described through an expres- 
sion with two parameters, simplified as E/E=(G[F/(EA)|°. 
Through correlation analyses, we connected the experimental- 
ly obtained coefficients, 1.e., the exponent a and amplitude (3, 
to the characterisation of rough surfaces. For a given material, 
the exponent a of this power law relation shows strong de- 
pendence on values of fractal dimension and RMS slope, 
while the coefficient @ is dominated by RMS roughness. 
These strong correlations can be used for establishing and 
validating predictive models for contact stiffness on a wide 
range of rough surfaces. 
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3.2 Paper 2: Numerical modelling of contact stiffness 


This second paper in this chapter presents a numerical approach for calculating the contact stiffness of 
fractal rough surfaces under compression. The numerical results are compared with experimental data, 
showing a quantitative agreement with measurements for a range of different rough surfaces. For 
diverse surface structures flattened by a rigid flat, a unified power-law function between the normal 
contact stiffness and the applied normal load can be found both experimentally and numerically. The 
exponent of the obtained power law relation increases with values of fractal dimension, while the 
magnitude decreases exponentially with the relative roughness amplitude. The numerical solutions of 
the proposed method are in good agreement with experimental results over a wide range of normal 
compression. The applicability of the presented method is further clarified in consideration of 
complex deformation of the compressed asperities under high stresses. Moreover, a parametric 
analysis is conducted to establish a correlation between contact stiffness and surface roughness 
descriptors. We found that the characterization method applied to determine the fractal dimension of a 
rough surface can affect somewhat the predicted contact mechanics. This study provides an easily 
implemented and computationally efficient method to connect mechanical behaviour with multi-scale 
surface structures, which can be utilized in design and optimization of engineering applications 
involving rough contacts, such as granular materials. The findings presented in Chapter 3 will be 
further linked with the electrical contact at rough interface presented in Chapter 4, in order to have a 
comprehensive understanding of the rough contact. The proposed numerical approach for rough 
contact along with the benchmark experiment presented in the first and second papers provides basis 
for modelling responses of granular materials under various stress states. 


This paper has been published in Jnternational Journal of Mechanical Science in 2017. I was the 
primary researcher and author of this paper, being supervised by Dr. Yixiang Gan and Dr. Dorian 
Hanaor. 
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In this paper, we study the contact stiffness of a fractal rough surface compressed by a rigid flat plane. A numerical 
model based on the analysis of flat punch indentation is proposed for simulated hierarchical surfaces, which 
are generated using statistical and fractal descriptors collected by surface profilometry. The contact stiffness of 
surfaces under increasing normal load is determined on the basis of the total truncated area at varying heights. 
The results are compared with experimental data from nanoindentation on four types of treated rough surfaces, 
showing good agreement with experimental observations below a certain truncation depth. Furthermore, the 
limits of the model’s validity are discussed by focusing on surface geometries and deformation of contacting 
asperities. With this proposed truncation method, we present a parametric analysis to establish a correlation 
between contact stiffness and surface roughness descriptors. The contact stiffness shows a unified power-law 
scaling with respect to the applied load over a wide range for simulated surfaces with distinct sets of roughness 
descriptors. The exponent of the power-law relationship is found to correlate positively to the fractal dimension 
while its amplitude is inversely correlated to the surface roughness amplitude. This study provides an easily 
implemented and computationally efficient method to connect mechanical behaviour with multi-scale surface 
structure, which can be utilized in design and optimization of engineering applications involving rough contacts. 


Keywords: 
Contact stiffness 
Rough surfaces 
Fractal dimension 
Nanoindentation 


1. Introduction 


The morphology of surfaces plays a determining role in interfacial 
phenomena including friction, adhesion, sealing, lubrication and ther- 
mal and electrical conductance [1,2]. When two solids with rough sur- 
faces are squeezed together, the area of true contact consists of numer- 
ous contact patches of various sizes and generally comprises only a small 
fraction of the nominal contact area. Over the years, numerous diver- 
gent approaches have been developed to explore load dependent con- 
tact behaviour at rough surfaces, yet the scientific community remains 
divided regarding the reliability of their predictions [1,3—6]. The vari- 
ous asperity-based approaches to contact models reported in the litera- 
ture can be categorized as: (1) multi-asperity contact models (i.e., sta- 
tistical models), in which the heights and / or curvatures of asperities 
follow given statistical distributions [7-10]; and (2) surface fractality 
models, including multi-scale models [11-17], the boundary element 
method (BEM) [6] and Persson’s theory [18, 19]. Further to this binary 
classification, the comparison between various reported models reveals 
discrepancies with noticeable differences between each other and with 
experimental results [17,19-21]. 

Within the above-mentioned approaches to surface structure simula- 
tion, it is also necessary to consider the contact behaviour of an individ- 
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ual asperity ranging from elastic, through elasto-plastic, to fully plastic 
deformation [8,22—29]. For purely elastic or plastic contact, the classic 
Hertzian model [30] and fully plastic model [31] can be applied, respec- 
tively. However, for elastoplastic regimes, individual asperity models 
may yield different contact responses, dominated by different deforma- 
tion mechanisms [32,33]. Particularly, a self-consistent analysis was put 
forward by Storakers, et al. [34] considering the general visco-elasto- 
plastic material indented by a spherical object. Additionally, shoulder- 
to-shoulder contact models for misaligned asperities were introduced to 
include oblique contact between pair asperities [15,35,36]. 

The contact behaviour of individual asperities can be combined 
to shed light on overall system behaviour by considering statistical 
and/or fractal approaches to describe surface morphologies. Pioneered 
by Greenwood and Williamson [7], multi-asperity contact models are 
based on the statistical height distribution (Gaussian or non-Gaussian), 
while assuming that the deformation of a given asperity is not influ- 
enced by that of neighbouring asperities. Within this framework, various 
implementations, considering different asperity geometries, have been 
developed over the past decades for the analysis of individual asperity 
deformation. In practice, statistical parameters for characterizing sur- 
face topography, such as variance of heights, slope, curvature, etc., have 
been used in this class of contact models. But these implementations as- 
sume features at a given narrow range of scales and thus depend on the 
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resolution of the surface measuring apparatus and sample length [13]. 
However, most natural surfaces exhibit features across a wide range of 
length scales [37], involving diverse morphologies, which bring about 
complexities in the modelling of interfacial properties [3,38]. In addi- 
tion, these statistical models were constructed by assuming that micro- 
contact forces arise principally due to deformation of asperities and are 
calculated considering only the base wavelength or a certain range of 
wavelengths, characterized by the roll-off and cut-off wavevectors of 
the power spectra of rough surfaces [4,12,17,19,39]. This assumption 
neglects the contribution of fine geometries at small length scales to the 
contact properties. In the fractal-theory related models, the finer sur- 
face features have been included in terms of overall contact area using 
an integration process [11,12,28,40]. 

An early appreciation for the significance of the multi-scale nature of 
surfaces was demonstrated by Archard [41]. In the contact theory using 
scale-independent parameters, Majumdar and Bhushan [42] utilized the 
fractal theory of Mandelbrot [43] to describe the distribution of contact 
areas between two rough surfaces. Yan and Komvopoulos [12] extended 
this theory to contact problems of three-dimensional rough surfaces, re- 
vealing the variations of the contact force and real contact area during 
quasi-static loading. Ciavarella, et al. [44] developed a two-dimensional 
fractal-based model for sinusoidal elastoplastic surfaces, and expanded 
the analysis to provide contact stiffness and interfacial resistance. Fur- 
ther relevant work was conducted by Jackson and Streator [13] using 
three-dimensional sinusoidal surfaces, considering the frequency spec- 
trum of the surfaces. Pohrt and Popov [6] deduced an empirical con- 
tact stiffness model by means of the boundary element method (BEM). 
The validity of the method of reduction of dimensionality, by which 
three-dimensional contact problems are mapped onto one-dimensional 
elastic contacts, has been investigated for non-adhesive contact of any 
axisymmetric bodies [45-47]. However, the fractal dimension has been 
shown to change due to an applied load [38,48,49] and thus an as- 
sumption of invariant fractal dimension in fractal-theory related mod- 
els may only be reasonable within a certain range of loading. Moreover, 
various methods [20,50,51] can be applied to quantify the fractality of 
a rough surface [19,20,52,53] possibly resulting in different values of 
the obtained fractal dimension, thus the adaptation of these methods 
merits further study to provide a consistent result in contact mechanics 
models. 

The true area of contact formed between two surfaces is of prime 
interest and is governed by morphology, material properties, and load- 
ing conditions. The contact properties of rough surfaces subjected to an 
applied normal load can be, at some level, interpreted by considering 
the true interfacial contact area along with surface descriptors [4,46]. 
However, determining the real contact area between contacting bod- 
ies through experimental or numerical analyses remains challenging. 
Factors that influence true contact area include the asperity height, cur- 
vature, the Poisson’s ratio of the material, strength and hardness and 
the existence of superimposed smaller asperities with similar or diver- 
gent scale-dependent properties [54]. With time, the sliding induced 
by expanding contacting spots [55] further affects contact morphology. 
Moreover, most natural surfaces exhibit features across a wide range of 
length scales the extent of which depends on simulation or measurement 
resolution [56,57]. 

Despite the fact that the contact area tends to present scale- 
dependent properties with contacting asperities in the status of incom- 
plete contact (with non-contacting zones surrounded by the contact- 
ing ones), the truncated areas at various depths for a given resolution 
can be employed to represent the real contact area. This was first pro- 
posed by Abbot and Firestone [31] to describe a wear process rather 
than indentation or flattening. Along this line, proposed truncation mod- 
els [58,59] have assumed that the contact area of an asperity pressed 
against a rigid flat can be approximately calculated by mathematically 
truncating the asperity tip. This is a reasonable assumption for small to 
medium interferences as within this range asperities tend to plastically 
deform due to the small radius of curvature. The average pressure be- 
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tween an asperity and a flat punch can simply be assumed to be equal 
to hardness, or can be related to material yield strength [7,26]. How- 
ever, Jackson and Green [22] showed that this simplification results in 
an inverse hardening process in which the hardness actually decreases 
with increasing interference. 

In addition to contact area the present work focuses on interfa- 
cial contact stiffness. An understanding of contact stiffness is impor- 
tant in contact mechanics, which plays a central role in governing the 
stress-dependent electrical and thermal transport between two contact- 
ing solids [6,60,61]. In the past decade, the relationship between surface 
structure and interfacial stiffness has been intensively studied numeri- 
cally and experimentally. Numerical analyses, using methods of molecu- 
lar dynamics, finite element analysis, etc., generally confirm linear pro- 
portionality between normal force and contact stiffness [39,62,63], as 
supported by the Greenwood-—Williamson model [7] and Persson’s the- 
ory [19,64]. However, other studies have reported that, for small to 
medium loads, the logarithm of stiffness exhibits close proportionality to 
the logarithm of the applied normal force [6,45,52]. In other words, the 
contact stiffness, k, is a power function of the normal force, Fy, ask « Fy, 
(with a<1), which differs from the work mentioned above. Numerous 
experimental studies have been carried out on rough surfaces to ascer- 
tain the relationship between surface structure and interfacial behaviour 
under load using diverse experimental approaches and materials. Jiang, 
et al. [52] measured the normal contact stiffness of cast iron specimens 
produced using different machining methods. Wang, et al. [19] mea- 
sured the contact stiffness of a rubber block squeezed against different 
concrete and asphalt road surfaces. Buczkowski, et al. [17] compared 
the normal contact stiffness determined using ultrasonic measurements 
with the fractal model based on Weierstrass—Mandelbrot function. Zhai, 
et al. [20] recently evaluated the contact stiffness at aluminium surfaces 
by nanoindentation tests utilizing different sized flat tips to achieve a 
wide range of applied stress levels. These experimental studies support 
the power law relationship between the contact stiffness and the applied 
load for certain stress ranges. 

The main purpose of the paper is to propose a comprehensive contact 
analysis method for a three-dimensional rough surface compressed by a 
rigid flat. The proposed truncation method is applied to simulated fractal 
surface structures characterized using various statistical and fractal de- 
scriptors, to interpret the variation of contact stiffness under increasing 
normal loading using experimental results for reference. Iteration proce- 
dures employed in simulating fractal rough surfaces ensure a description 
with identical parameters which are critical in determining the normal 
contact stiffness. The applicability and repeatability of the presented 
method is discussed based on geometrical features and analyses of the 
studied rough surfaces. This study provides an easily-incorporated and 
highly-effective numerical method for predicting contact stiffness un- 
der conditions of small to medium loads. Following validation, a para- 
metric analysis is conducted for simulated surfaces varying in surface 
topologies, allowing the obtained contact stiffness to be related to sur- 
face roughness descriptors. Finally, correlations between normal contact 
stiffness and roughness descriptors have been established and discussed 
with respect to fractal dimension values obtained using different meth- 
ods. 


2. Theoretical framework 


The deformation mechanics of contacting surface asperities remains 
a topic of significant debate in the research community. Under most 
applied conditions the true contact area between rough surfaces involves 
only a small fraction of the nominal contact area, and for this reason the 
total real contact area can be considered to approximately equal to the 
truncation area in the present method. We consider a simple extension of 
the contact analysis for indentation by a flat punch [65], which reveals 
that there is a relation between contact stiffness, contact area, and elastic 
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modulus, given by 


EEE) =P : 

r ©? i 5 VA 
where K is the contact stiffness in units of N/m, A is the apparent con- 
tact area and f is a geometrical constant, taken as unity for a flat punch. 
E, denotes the reduced Young’s modulus including contributions from 
the both the compressed rough surface and opposing flat surface, given 
by 1/E, = (1 — v2)/E, + (1 — v?)/E,, with the subscript i indicating the 
opposing surface, the subscript c indicating the properties of the com- 
pressed surface and v being the Poisson’s ratio. The contact force can 
be obtained through integration of the contact stiffness with respect to 
the indentation increment. Eq. (1) is a fundamental equation for assess- 
ing the elastic properties in indentation tests, and has been shown to be 
equally applicable in cases of elastic-plastic contact. If the elastic mod- 
ulus is known, the relation between contact area and contact stiffness 
can thus be obtained [66,67]. 


(1) 


2.1. Validation for single asperity under compression 


An asperity under compression can behave elastically or plasti- 
cally. For describing purely elastic contact of rough surfaces, the classic 
Hertzian model can be applied as suggested by previous contact mod- 
els [7,26,33,42]. For a sphere under compression by a rigid plane, it 
is assumed that the radius of a contacting asperity, R, is very close 
to that of the undeformed case. In such cases, the relationship be- 
tween R and the truncation microcontact radius, r’, can be obtained by 
R* =(R—d)*? + (r’)?, where d is compression depth. Since R is typically 
orders of magnitude greater than the d, the relationship can be reduced 
to (7’)? = 2Rd [12,27]. In Hertzian contact, the radius of the real contact 
area is given by r = \/ Rd. Thus for a circular microcontact in the fully 
elastic regime the relationship between the truncation area (a’) and the 
real elastic contact area (a;) can be expressed as: a’ = 2a,. 

For describing plastic deformation, Abbot and Firestone [31] devel- 
oped the most widely used model for a fully plastic contact, known as the 
surface microgeometry model. This model assumed that the deformation 
of a rough surface against a rigid flat plane is equivalent to the trunca- 
tion of the undeformed rough surface at its intersection with the plane. 
As a result, the real area of the contact can be approximated simply as 
the geometrical intersection of the flat with the original undeformed 
surface profile, and the contact pressure is the plastic flow pressure. In 
addition, McFarlane and Tabor [68] showed that both the normal and 
tangential stresses contribute to the deformation of junctions formed at 
the interface of contacting bodies under compression. For cases of elastic 
contact, the tangential stress constrains the expansion of the compressed 
asperity through static friction. As the contacting region between an as- 
perity and the rigid flat is already in the plastic state of stress under 
purely normal loading, tangential loading, however small, may cause 
further yielding. Equilibrium can be maintained if the area of a contact- 
ing spot grows [68]. Therefore, for a fully plastically microcontact, the 
relation between the truncation area and the real plastic contact area 
(ap) is: a’ < Gp. 

The above discussion supports that the truncation area, a’, of asingle 
compressed asperity at any given interference lies in the range of ag < a’ 
< dp, the upper and lower limits of which are determined by the contact 
area for fully elastic and plastic regimes, shown in Fig. 1. For a single as- 
perity of a fractal rough surface, exhibiting a hierarchical structure, we 
also assume that the contact area of an individual asperity for a given 
surface interference will follow the relation: ag <a’<ap. Furthermore, 
the truncation method could provide contact properties between elastic 
and plastic behaviours. This assumption is further rationalized by com- 
paring the dependence of contact area on contact force obtained from 
the concise self-consistent model for contact of inelastic materials, based 
on the indentation theory and von Mises isotropic flow theory [9]. 
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For a visco-elasto-plastic contact problem, the constitutive behaviour 
can be simplified as o = ope”2é"8, where og is a material yield param- 
eter and mg and ng are hardening and creep exponents, respectively. 
For a time-independent perfectly plastic material, mg =ng =0, while in 
the linearly elastic case, where mg = 1 and ng =0. A semi-analytical so- 
lution for the contact problem can be obtained for spherical bodies [9, 
34]. In this paper, only a quasi-static situation is studied without con- 
sidering creep, i.e., Ng =O. Here, we consider a spherical asperity under 
compression by a rigid plane using the proposed truncation method and 
benchmark this special case to other solutions existing in the domain of 
contact mechanics, in Fig. 2. The maximum interference is 100 times 
the critical interference, w, =(zkH/2E,)2R where k=0.454 + 0.41v is 
the hardness coefficient of the asperity related to its Poisson’s ratio, v, 
and, H is the hardness of the asperity relative to its yield strength. This 
critical interference value marks the inception of the elastoplastic defor- 
mation [13,26,27]. The elasto-plastic solutions for materials exhibiting 
different hardening exponents (mg =0.0, 0.10, 0.20, 0.30, 0.40, 0.50 
and 1.0) are shown in Fig. 2. The obtained numerical results from the 
presented truncation method (1000 truncation depth increments) lie be- 
tween fully plastic and elastic regimes. Thus, to a certain extent, these 
results demonstrate the applicability of the truncation method for the 
analysis of elasto-plastic contact. 


2.2. The truncation method for rough surfaces 


In this work an approximate model is presented for estimating the 
contact stiffness of a rough surface compressed by a rigid plane with the 
following assumptions: 


(1) The studied fractal rough surface is flattened by the rigid plane, 
and contact ‘islands’ grow in terms of size and number through 
successive truncations parallel to the mean surface plane; 

(2) The truncated surface features on the rough surfaces flow plasti- 
cally into the valleys of non-contacting regions; 

(3) Contacting asperities do not interact with each other; 

(4) Deformation is confined within the interface region rather than 
the bulk region. 


In this method, a rough surface is levelled through a mountain-top re- 
moval approach in order to obtain truncation areas. For a single asperity 
under compression, it is a reasonable assumption to use the total trun- 
cation area, A! (Ap < A’ = Y=" a, = A, < Af, where A, is the true con- 
tact area), at a given height to approximate the real contact area for the 
corresponding surface interference, w. The contact stiffness is extracted 
based on the analysis of indentation tests by a flat punch [65]. Under 
increasing compression, smaller microcontacts expand, merging to form 
larger ones. Consequently, the value of contact stiffness, E. depends pri- 
marily on the true contact area and approaches the elastic modulus of 
the bulk material, E. The contact stiffness (k, with the unit of N/m) can 
be estimated by the following expression using the true contact area 
(A.), based on Eq. (1): 


k=pLEVA, (2) 
ve 

where E’ is the (constant) reduced elastic modulus calculated for the 
bulk elastic properties of the tested material (Young’s modulus, F and 
Poisson’s ratio, v) and indenter: E! = [(1 — v”)/E + (1 —0?)/E,]"!, and 
f’ is a geometric factor of the order of unity, equalling to 1. By writing 
Eq. (2), we assume that the effect of surface roughness on the measured 
incremental stiffness can be described by considering the true contact 
area A, (rather than the projected contact area A) and bulk material 
properties in Eq. (1). By comparing Eqs. (1) and (2), and considering 
that E; > E > E., we obtain the following scaling relation: 


E../E* « piv A./A. (3) 


Considering dF=k(w)d@, we can then obtain the loading force F 
through integration of contact stiffness, k(w), during the compression 
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Fig. 1. Schematic of a deformed asperity under compression by a rigid plane. The areas of a, and dp, show the real contact areas for fully elastic and plastic regime, respectively, while 


a’ represents the truncation area for a given surface interference of w. 
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Fig. 2. Comparison of obtained contact responses of a single spherical asperity under 
compression of a rigid flat between the truncation model and the inelastic flattening model 
[34]: The normalised contact area versus normalised load is shown with A equalling zR?, 
E being the Young’s modulus of the material. 


with penetration increments da. Note that we are here focusing on the 
hierarchical structure of the surfaces by applying normalization proce- 
dures. 


3. Surface characterization and simulation 


To gain insights into the validity of the presently developed trun- 
cation method for rough surfaces, we experimentally and numerically 
evaluated the contact stiffness for round disk samples made of alu- 
minium alloy with surface treatments applied in order to produce dis- 
tinct surface structures. 


3.1. Surface fabrication and profilometry 


To produce different rough surfaces, three surface treatment meth- 
ods were applied: (1) polishing, (2) surface mechanical attrition treat- 
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ment (SMAT), and (3) abrasive blasting using various sized glass beads. 
More details of the surface treatment and topographical characterisation 
can be found in [20]. Fig. 3 shows typical scanning electron microscope 
(SEM) images, 3D surface topography and corresponding simulated sur- 
faces. Over 10 digitized scans, the treated sample surfaces were char- 
acterised through descriptors of the peak-valley height R,, root mean 
square roughness R,,,,, roll-off wavelength W,, and fractal dimension 
D characterized through different methods, as shown in Table 1. These 
roughness descriptors will be used later to generate simulated rough 
surfaces that have similar features for evaluating the proposed trunca- 
tion method. Using the fractal dimension we are able to estimate sur- 
face structures at scales below the equipment resolution, assuming self- 
similar scaling in terms of geometrical features [37]. This estimation 
approach further improves the efficiency of the present method as utilis- 
ing surface scans at sub-micron scale resolutions, such as those obtained 
through AFM, would be computationally intensive and would provide 
limited information from higher scale features. 


3.2. Simulated rough surfaces 


The scale-invariant parameter, i.e., fractal dimension, provides a 
means of describing realistic multiscale roughness. Previous studies 
[12,52,55] have shown that a fractal rough surface can be determin- 
istically simulated by the Weierstrass—Mandelbrot fractal function [43], 
which can be written as 


(D,,-3) 1/2 M Max 
z(x,y) = wr(<) (42) 2 2 ( (Par 3) 


1/2 
x< cos cos mal (x° = y) cos (tan~( *) 7 ) — 
a, — a = a 
ee W, x M me 
(4) 


The parameter y determines the density of frequencies used to con- 
struct the surface and, in similarity to previous work [33], is set as 1.5 
based on considerations of the surface flatness and frequency distribu- 
tion density. G is a height scaling parameter independent of frequency, 
termed topothesy [48]. The parameter W.. is the roll-off wavelength, 
which can be obtained from the power spectrum, determining the ba- 
sic wavelength of highest scale features. In this paper, values of roll-off 
wavelength, W,., of the treated surfaces are obtained by estimating the 
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Fig. 3. SEM images, typical digitized profilometry scans (1024 x 1024 pixels over an area of 1 x 1mm?) of aluminium samples subjected to different surface treatments and their 
corresponding simulated surfaces (10, 240 x 10, 240 pixels) presented in row 1, 2 and 3, respectively: (a) polished; (b) SMAT with 2 mm-sized steel balls; (c) blasted with 300 um-sized 


glass beads; (c) blasted with 50 um-sized glass beads. 


Table 1 


Surface characterisation for digitized scans and simulated surfaces (indicated with a suffix ‘s’ at the end of corresponding sample names). 


Fractal dimension 


Dri 


Fractal dimension Fractal dimension Fractal dimension 
D box Dys Dsp 


2.093 +0,0618 


Surface Amplitude Root mean square Roll-off 
roughness R; or roughness wavelength 
R,’ /um Rrms /¥M W,/pm 
Polished 3.20341,939 0.057 + 0.005 > 1024 
Ps 0.263 + 0.032 0.053 + 0.004 > 1024 
SMAT 22.184 +5.201 2.730 0.255 ~120 
SMATs 16.2104 1.569 2.762 + 0.053 120 
GB300 42.376 + 9.238 4.179 + 0.194 ~ 60 
GB300s 36.865 + 1.731 4.102 + 0.034 60 
GB50 32.239 + 7.843 2.970 + 0.276 ~25 
GB50s 2.935 + 0.031 25 


2.091 + 0.0344 
2.228 + 0.0199 
2.229 + 0.0232 
2.991 + 0.0217 
2.555 + 0.0425 
2.626 + 0.0174 
2.625 + 0.0328 


2.024 + 0.0410 
2.025 + 0.0313 
2.15620,0131 
2.161 + 0.0297 
2.424 + 0.0257 
2.417 + 0.0316 
2.351 + 0.0363 
2.391 + 0.0395 


2.103: +0,0628 
2.103 + 0.0562 
2.422 + 0.0436 
2.417 + 0.0487 
2.811 + 0.0769 
2.813:4-0,0521 
2.824 + 0.0767 
2.825 + 0.0876 


2.317 + 0.0788 
2.315 + 0.1196 
2.412 + 0.0964 
2.421 + 0.1222 
2.284 + 0.0931 
2.282 + 0.0876 
2.279 + 0.1054 
2.284 + 0.0933 


26.239 4 1.157 


average distances between asperities, based on surface structures mea- 
sured experimentally. The frequency index n,,,, assumes finite values 
corresponding to the cut-off wavelength, W,, representing the small- 
est distance between two adjacent pixels [12,55]. The parameter M de- 
notes the number of superposed ridges used to construct the surface. 
A random number generator is used to uniformly distribute the val- 
ues of random phase a,,,,. Here, isotropic fractal surfaces, consisting 
of 10,240 x 10,240 pixels, with M= 10, were generated to simulate real 
surfaces based on the surface descriptors obtained from digitized scans 
with the projective area being 1 mm x 1mm, shown in Fig. 3. 

In the fractal analysis of surface asperities, the fractal dimension and 
the topothesy as are commonly considered as two invariants. However, 
these two roughness parameters have been shown to vary with load 
[48]. In this paper, a scaling procedure based on amplitude roughness, 
R,, was employed to govern the vertical properties of simulated surfaces 
instead of topothesy, G (set as 1). The values of all pixels, representing 
the height, were scaled according to the height-width ratio (R,/L, where 
L is the side length of the surface) of the simulated surface, in order 
to yield the target peak-valley height. The fractal dimension of rough 
surfaces (for both digitized scans and simulated surfaces) was then es- 
timated using methods of (1) triangulation [20,69], (2) box-counting 
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[70,71], (3) vertical sections [51], and (4) power spectrum analysis 
(Mitchell and Bonnell, 1990; Williams and Beebe Jr, 1993; Van Put et 
al., 1994; Mannelquist et al., 1998), with the obtained fractal dimension 
marked with D,,;, Dyox, D,;, and D,,, respectively. The four methods em- 
ployed are widely used in the calculation of fractal dimension [61,72- 
74], though the resulting fractal dimension may be not identical. As is 
shown schematically in Fig. 4, the procedure for obtaining the input 
fractal dimension D,, for simulated surfaces involves assuming an ini- 
tial value of D,,, (e.g., the value obtained from the real surface scan). A 
surface (1024 x 1024 pixels over an area of 1 x 1 mm7) is generated with 
inputting other required parameters, i.e., the roll-off wavelength W,. and 
amplitude roughness R,. Subsequently, the calculated value of fractal 
dimension, D,,,, of the generated surface is compared with the value 
obtained from the real surface scan and D,,, is adjusted accordingly. 
This procedure is repeated until the calculated fractal dimension, D,,;, 
approaches that determined from the profilimetry of real surfaces using 
one chosen method (e.g., error < 5%). Even though the triangulation, 
box-counting method, vertical sections approaches and power spectrum 
analysis result in different values of fractal dimension, a larger input 
fractal dimension consistently leads to a larger evaluated D,,,, with all 
four approaches. 
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Fig. 4. Flow chart of iterative procedures for surface simulation. 


As amplitude roughness R, from digitized scans is sensitive to mea- 
surement resolution and noise, R,,,, is employed to yield vertical fea- 
tures approximating real surfaces. Similar to the fractal dimension, an 
iteration procedure is employed to obtain a proper R,’ value with which 
a simulated surface can replicate the vertical geometrical features of 
the digitized scans. Subsequently, the obtained D,,,’ and R,’ from the it- 
eration procedures, along with the basic wavelength, W,, are used to 
generate rough surfaces at higher resolutions. The sequence of the two 
iteration procedures, for D,,,’ and R,’ should not be reversed, as the vari- 
ation of input fractal dimension can influence the roughness amplitude, 
but not vice versa. The simulation of surfaces using the proposed itera- 
tion procedures are presented in Fig. 4. 

The amplitude roughness, R,, and the RMS roughness, R,,,,, have 
been shown to be scale-dependent, which can be determined by the 
roll-off and cut-off wavelength [11]. The values of the two parameters 
used in the iteration procedure, shown in Fig. 4 and Table 1, are those 
calculated at the scanning resolution. 


3.3. Influences of resolution 


As the true contact area between real surfaces is difficult to defini- 
tively determine, it is important to investigate the influences of the sur- 
face resolution, as well as step sizes of truncating increments. The to- 
tal contact area can be determined by summarising the contribution of 
contacting asperities. For a given truncation, the total truncated area, is 
decreasing as the resolution increases [12,43]. 

The contact stiffness obtained using the proposed method was stud- 
ied using simulated surfaces of varied resolution. Results obtained from 
scans with resolution ranging from 1024 x 1024 to 32,760 x 32,760 pix- 
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els over an area of 1x 1mm? are compared in Fig. 5(a). The effects of 
truncation increment size are shown in Fig. 5(b). For both scenarios, the 
maximum truncation depths are up to half of the roughness amplitude, 
which is roughly the mean plane of the rough surface. The shaded error 
bars are obtained from ten simulations for each resolution and trunca- 
tion step size. For relatively low-resolution surfaces, a clear resolution- 
dependence can be observed, while this trend tends to diminish as res- 
olution increases. Similarly, simulations using 200 or more increments, 
up to 5000, appear to have small differences. Therefore, a resolution 
higher than 8192 x 8192 with more than 200 increment steps appears 
sufficient for the truncation method to provide results yielding negligi- 
ble resolution sensitivity and numerical convergence. 

With an increasing resolution, the contact area, and therefore, the 
contact stiffness approaches the asymptotic limit. This is consistent 
with the other fractal-based models [3,4,13,14,18], in which scale- 
independent contact area and contact stiffness can be achieved. Even 
though a converging trend has been observed based on the proposed 
framework, the question remains on quantitatively estimating the con- 
tact area, interfacial void, contact resistance, and any other quantity 
which depends on extremely fine details. The ongoing fractality shown 
by the hierarchical surface structures may effectively become less impor- 
tant for interfacial properties on a critical length scale marking the de- 
fined mechanical interaction, electron transport, heat transfer, etc. This 
can be generally supported by experimental research with observed val- 
ues of interfacial parameters, such as contact stiffness, electrical contact 
resistance, and friction etc. [3,20,52,61,62,75-77]. This convergence 
trend will be discussed in details in our following work. In this paper, the 
horizontal resolution of the generated surface is set as around 100 nm, 
i.e., 10,240 x 10,240 pixels over an area of 1X1mm7, and 200 incre- 
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Fig. 5. Sensitivity studies on spatial resolution and increment steps on the calculated contact stiffness using the truncation method: (a) pixel level; (b) number of increment steps. The 
inputs of the compressed surfaces are as follows: D,,= 2.5, R,x= 10,000 nm and W,=100 um. 


ment steps are considered in the integration process for estimating the 
contact force. 


4. Applications and discussion 
4.1. Nanoindentation on roughened surfaces 


The normal contact stiffness of three types of treated surfaces was 
assessed using nanoindentation (Agilent G200) with three flat inden- 
ter tips of different diameters of 54.1 um, 108.7 um and 502.6 um (FLT- 
D050, FLT-D100, and FLT-D500, respectively, SYNTON-MDP). In order 
to evaluate only the elastic responses, partial unloading tests were suc- 
cessively performed at ten intervals by decreasing the applied load by 
10% each time. The loading level of each successive unloading stage 
is twice that of the preceding one, with a maximum load of 500mN 
reached during the final unloading step [20,61]. The obtained contact 
stiffness, in units of N/m was further transformed to surface stiffness in 
units of Pa through E, = (4/7K) if (2A) where A is the projected area 
of the flat indentation tip. The experimental results show that the re- 
lation between the contact modulus and applied loading force can be 
described as E./E=f[F/(EA)]%. Correlations between the surface pa- 
rameters and, a and f the evaluated are detailed in [20]. From the 
present experiments, the exponent a of the obtained power law relation 
shows strong dependence on surface fractality, while the roughness am- 
plitude is found to govern the magnitude of parameter f£. The maximal 
indentation depth reached is around 5 um, which is comparable with the 
roughness amplitude. The experimentally observed low-level compres- 
sion depths suggest that the contact stiffness over a wide stress range 
is strongly affected by finer surface details, which are superimposed on 
the basic waves. These finer features may not be included in the study of 
contact mechanics governed by the surface structure [3—5,12]. Specif- 
ically, the considered finest feature of the roughness can be different 
in various engineering applications. The electrical contact resistance 
[61,76] and adhesion [78] have found to be sensitive to fine surface 
structures at smaller length scales. However, the influences of surface 
structures on interfacial thermal conduction can be well described using 
RMS surface roughness, which depends primarily on basic wavelength 
[79]. 
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4.2. Comparison with experimentally measured contact stiffness 


Contact analyses using the proposed truncation method are pre- 
sented in Fig. 6. For SMAT and GB surfaces, the numerical results 
obtained from simulated surfaces are in good agreement with exper- 
imental results, particularly for low applied pressures, correspond- 
ing to small truncation depths, in the range 0 to 10.0 ~ 15.0%R,’, 
30.0 ~ 35.0% R,’, 32.5 ~ 37.5% R,’, and 25.0 ~ 30.0% R,’ for simulated 
surfaces, Ps, SMATs, GB300s, and GB50Os, respectively. 

It can be found in Fig. 6 that the predicted results for polished sam- 
ples show the most significant divergence between numerical and exper- 
imental results in particular for interferences greater than 15% R,’. This 
is likely the result of several factors, notably the misalignment of the in- 
denter with tested surfaces and roughness features that may be present 
on the flat diamond tip (measured to have an RMS roughness of around 
0.05 um) as well as adhesion, boundary effects and surface hydration, 
all of which yield a more pronounced effect for polished surfaces. Even 
though differences between the numerical analyses and experimental re- 
sults can be noticed, by comparing the numerical results obtained from 
the mentioned models with the experimental results, we can systemati- 
cally deduce the following conclusions: (1) power-law relation between 
the contact stiffness of a fractal rough surface and the normal load is 
found; (2) the roughness amplitude represented by R, and R,,,, plays 
a considerable role in determining the amplitude of the function; (3) 
the slope of the obtained power functions of the contact stiffness on 
applied normal force are in the range [1/3, 1), depending mainly on 
the fractal property of the surface. Comparative analyses of the exper- 
imental results with existing models can be found in [16,17,20,52,61]. 
It worth noting that linear relationship can be obtained within a certain 
range of normal compression, as shown in [17,19,64,80]. More impor- 
tantly, in this paper, power-law relationships over a wide range of nor- 
mal compression load have been experimentally observed and numeri- 
cally demonstrated using our proposed method, covering the small-to- 
medium load. 


4.3. Applicable truncation depth 


As shown in the previous section, the proposed truncation method is 
capable of efficiently extracting contact stiffness on the basis of surface 
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Fig. 6. Comparison of the experimental results and the contact stiffness obtained using the truncation method using D,, . Normalised stiffness, E-/E, and normalised applied load, F/(EA), 
are used, with E being the Young’s modulus of the tested material, and A the apparent contact area. The numerical results based on the present method are shown using dashed lines 
with shaded error bars. Solid lines (representing the fitted power functions with the exponent values provided) with error bars are for experimental results. 





Fig. 7. Truncation sections for a typical scanned GB surface at varying heights. The truncation starts from the top of the surface with depths of 10%R,, 20%R,, 30%R,, 40%R,, and 


50%R,, and the estimated contact areas are coloured with red. 


descriptors for small to medium loads. This method provides an easily 
implemented numerical framework to relate contact behaviour to the 
geometrical description of studied surfaces, with tunable resolution and 
calculation accuracy, meeting the requirements for various engineering 
applications. However, as the load and the truncation depth increase, 
the applicability of the presented method merits further discussion, due 
to the complex deformation of the compressed asperities, interaction 
between neighbouring asperities, adhesion, and friction. The applicable 
interference range is discussed below by concentrating on the surface 
geometries and asperity interactions. 

In this model, for a given surface interference, the total contact area 
under compression is simply determined from numerical integration of 
the truncated areas of the rough surfaces, at the corresponding trunca- 
tion height. As shown in Fig. 6, the proposed numerical procedure is 
capable of predicting the contact stiffness until the truncation depth 
reaches a certain value. A series of truncation sections for a typical 
GB300s surface are shown in Fig. 7. As the surface interference in- 
creases, small contacting spots expand, giving rise to increased inter- 
actions between asperity contact regions. A uniform rise in the height 
of the non-contacting regions has been experimentally observed [81], 
showing interactions also at asperity bases. Increasing compression can 
also bring about a variation in asperity shape. As the surfaces approach, 
contacting spots continue to grow and merge with each other, with me- 
chanical behaviour transitioning gradually to that of a bulk material. 
These phenomena at large surface interferences result in limitations to 
the present truncation method beyond a certain level of deformation. 

In order to quantitatively define the applicable range of the trunca- 
tion method, we firstly consider the variation of island numbers with 
respect to depth, which can be employed to indicate the interaction be- 
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tween asperities. The numbers of contacting asperities at various trun- 
cation depths for the four types of studied surfaces over a projected 
area of 3x 3mm? with 30,720 x 30,720 pixels are shown in Fig. 8(a). 
This size is chosen to give a meaningful statistical representation of the 
islands. As the interference depth increases, the number of contacting 
spots reaches a maximum at depths larger than those maximum trun- 
cation depths. The range of truncation depths in which the growth rate 
of island numbers increases monotonically coincides with the applica- 
ble range of the present method, as shown in Table 2. Furthermore, we 
considered the evaluation of perimeters and areas of all levelled islands 
at varying truncating heights using the approach proposed by Mandel- 
brot (1984). The correlation law between perimeter, L;, and Area, A, 
for sections of truncated asperities, like islands surrounding by water, is 
given by: 

bp eA Pe? (5) 
Within a certain range, the correlation between perimeter and area fol- 
low a power law type behaviour for all the four types of studied surfaces. 
However, as the truncation depth further increases, the exponent expe- 
riences a downward trend as is shown in Fig. 8(b), departing from the 
power law behaviour shown in Eq. (5). This suggests a critical trunca- 
tion depth range, in which a truncated surface can be described by an 
invariant fractal dimension. In this paper, at any given truncation depth, 
the curve fitting between L; and A are conducted, with all obtained data 
for L; and A to extract the exponent value, and the critical depth is ob- 
tained as an obvious decrease for the exponents is observed, as listed in 
Table 2. The exponents of the obtained power law relations are all larger 
than 0.5, which is the typical value of the exponent obtained in verti- 
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Fig. 8. Surface geometry analyses: (a) normalised numbers of contacting spots at various truncation depths for the studied surfaces, where the critical truncation depths, corresponding 
to maximal growth rates of island numbers, are marked with dashed lines; (b) typical log-log plots based on slit island analysis for the four types of studied surfaces (four simulations 
have been carried out for each type of surfaces) with the solid lines having a respective the exponent value, shown in Eq. (5). 
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Fig. 9. Dependence of the normal contact stiffness on the normal pressure: (a) for constant fractal dimension, D,,, = 2.5 and different relative roughness 6; (b) constant relative roughness 


6=0.01 and different fractal dimension D,,;. 


Table 2 
Indicators for the applicable range of the proposed truncation method. 


Surface Maximum applicable depth for truncation 

method (% R;) 
Ps 10.0-15.0 33.0-33.5 
SMATs 30.0 -35.0 36.0-36.5 
GB300s 32.5-37.5 31.0-31.5 
GB50s 25.0-30.0 34.0-34.5 


cal section analysis for an asperity in regular shape, such as a sphere or 
cone. 

Despite the fact that the two indicators we choose here are based only 
on surface topology, these geometrical features are closely related to the 
mechanical responses under normal compression. Within the range of 
truncation depths suggested by the two indicators, contacting asperities 
tend to be sufficiently isolated and thus respond to the external load- 
ing individually. Moreover, in this range, an invariant fractal dimension 
can be applicable to describe the multi-scale nature of the surface. At 
a truncation depth beyond this range, interactions between neighbour- 
ing asperities will have a significant influence on the contact, causing a 


Depth for the maximum increasing rate of 
island number (% R;) 


313 


Maximum depth for surface exhbiting 
fractality (% R;) 


29.0-30.0 
32,.5-37.9 
30.0-35.0 
21 -o200 


change of asperity shape, non-ignorable adhesion, sliding at interfaces 
etc. Even though the applicable range may also be affected by the me- 
chanical properties of the material, including hardening, Poisson’s ra- 
tio, and scale-dependent yielding strength, the purpose of this paper is 
to present a numerical framework to consider the structure dependent 
contact behaviour. 


4.4. Parametric study using the truncation method 


In order to relate the surface structure, represented by descriptors 
including fractal dimension, roughness amplitude and wavelength, to 
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Fig. 10. (a—e) Correlation between the stiffness variation exponent a and values of fractal dimension characterized using different methods; (f) Correlation between the stiffness magnitude 
6 and values of roughness amplitude, 6. Horizontal and vertical error bars show the corresponding standard deviations from the surface measurement and curve fitting. 


contact stiffness, contact analyses are conducted for fractal rough sur- 
faces with various surface topographies. The contact analyses were per- 
formed under conditions of small surface interferences as discussed in 
the previous section, in a regime where the truncation method performs 
well. Two groups of surfaces (10,240 x 10,240 pixels over an area of 
1x1mm7) were generated and analysed by the present model, with 
the maximum surface interference being 0.1R,. One group of simu- 
lated surfaces was generated with a constant fractal dimension value 
of D;, = 2.50 and a varied roughness amplitude 6 (0.000625, 0.00125, 
0.0025, 0.005, 0.01, 0.02, 0.04, 0.08, and 0.16). For the second group, 
surfaces were set with the same relative roughness amplitude, described 
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as 6 =R,/W,=0.01 and exhibited varying fractality D,,, (2.10, 2.30, 2.50, 
2.70, and 2.90). In order to further test the repeatability of the present 
model, we analysed five simulated surfaces for each individual set of in- 
put roughness parameters. Fig. 9 shows that the predicted contact stiff- 
ness exhibits a power-law relation with the applied normal load, which 
can be described as E./E= f[F/(EA)]*. The parameters of exponent a 
and magnitude f are related only to surface structure. The exponent a 
increases with values of fractal dimension, D,,,, while the relative rough- 
ness amplitude, 6, indicating primarily the vertical scale of the rough 
surface, dominates f. The correlations for the stiffness parameters, i.e., 
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Table 3 


Fitting functions for the stiffness parameters, a and f, using D,, and 6. 


Fractal dimension a 

Din a = 0.209D,,, + 0.090 0.940 
Di a = 0.313D,,, — 0.186 0.948 
Dyox a = 0.394D,,,. — 0.347 0.966 
D,, a =0.328D,,, — 0.258 0.980 
Dsp a=0.229D,, + 0.071 0.741 
5 \ 


a and f, using 6 and D were obtained and shown in Fig. 10. The corre- 
sponding fitted functions are presented in Table 3. 

The plots shown in Fig. 10 examine the effects of using different 
methods for the evaluation of the fractal dimension of rough surfaces 
in contact [3,6,11-13,20,61], namely the triangulation method (D,,,), 
box-counting method (D,,,.), vertical sections method (D,,), and power 
spectrum analysis (D,,). The obtained fractal dimensions are expected 
to differ somewhat due to the use of different calculation methods. As 
shown in Fig. 10 and Table 3, using fractal dimension obtained from dif- 
ferent methods will produce different fitting functions in terms of a(D). 
Thus it is necessary to clarify how the fractal dimension is extracted and 
the applicable interference and/or load range of the chosen approach of 
contact mechanics. However, the regression factors for fitting functions, 
a(D), are greater than 0.94 in most cases, except the one using the power 
spectrum analysis, demonstrating a strong correlation with fractality. 


5. Conclusion 


In this paper, we propose an efficient numerical approach for calcu- 
lating the contact stiffness of fractal rough surfaces under compression 
by considering a series of truncation sections and discuss the applica- 
ble range of this method. Numerical predictions of contact stiffness are 
benchmarked with experimental data over a wide range of applied nor- 
mal loads, showing a quantitative agreement with measurements for a 
range of different rough surfaces. Parametric analyses were carried out 
for simulated surfaces with varying surface topologies, to establish cor- 
relations between the obtained contact stiffness and surface roughness 
descriptors, including the fractal dimension and roughness amplitude. 
Based on the results and discussion presented, the following conclusions 
can be drawn: 


(1) Fractal rough surfaces can be well described using three key pa- 
rameters including roll-off wavelength, amplitude roughness and 
fractal dimension. For diverse surface structures flattened by a 
rigid flat, a unified power law function between the contact stiff- 
ness and the contact load can be found both experimentally and 
numerically. The exponent of the obtained power law relation 
increases with values of fractal dimension, while the magnitude 
decreasing exponentially with the relative roughness amplitude. 
The numerical solutions of the proposed method for contact stiff- 
ness versus contact load are in good agreement with experimental 
results over a wide range of normal compression. The surface in- 
terference range over which the proposed truncation method is 
applicable, can be determined by monitoring the variation of the 
number of contact islands in truncated sections, and the evolu- 
tion of the total perimeter and area of these islands, as trunca- 
tion deepens. These geometrical features are closely related to 
the mechanical responses under normal compression indicating 
the interactions between asperities. Within the applicable range, 
an invariant fractal dimension can be employed to present the 
hierarchical properties of the surface structure. 

The parametric analyses presented in this study suggest that dif- 
ferent methods for characterizing the fractal dimension of a rough 
surface may lead to different predicted contact stiffness within 
the applicable interference for approaches based on fractal the- 
ory. This implies that the characterization method applied to de- 


(2) 


(3) 


Coefficient of correlation (for a) 
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B Coefficient of correlation (for B) 
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=0,03865~-""" 0.992 


termine this surface descriptor can affect somewhat the predicted 
contact mechanics. 
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4 Electrical-mechanical Behaviour at 
Rough Interfaces 


Electrical contact is a complicated phenomenon due to the multi-scale surface structures and many 
involved mechanisms of resistance modification, including dielectric breakdown, localized current- 
induced welding, percolation collective process, and arising chemical reactions. The understanding 
the involved multi-physics coupling between electrical, thermal and mechanical fields is of 
considerable importance in quantifying the single contact resistance at a rough interface and the 
overall contact macroscopic resistance throughout the granular media. 


At 2 


4.1 Paper 3: Electrical contact resistance at rough interfaces 


This paper presents an experimental investigation on the electrical contact resistance between rough 
surfaces in close contact under various compressive stresses. We considered isotropically roughened 
aluminium disks with upper and lower surfaces modified through polishing and sand blasting using 
different sized glass beads. Fractal geometry and roughness descriptors, including root mean square 
values of roughness and slope, were used to describe the topography of sample surfaces, based on the 
digitized profiles obtained from interferometry-based profilometry. The electrical contact resistances 
at the interfaces were measured through various scenarios realised by controlling experimental 
parameters including testing time, testing current and mechanics loading. The experimental results 
show that the measured resistance depends closely on the measurement time, testing current, surface 
topology, and mechanical loading. For properly chosen measurement time and testing current, we 
demonstrated that a power-law relation exists between the measured electrical resistance and normal 
stress across a certain range, with the absolute values of exponents rising with the fractal dimension of 
the surfaces. The work will be of interest to a general community in contact mechanics, as well as in 
material engineering. The observed evolution of electrical contact resistance at various loading 
conditions helps to understand the electrical transport in granular materials. 


This paper has been published in ASCE Journal of Engineering Mechanics in 2015. I was the primary 
researcher and author of this paper, being supervised by Dr. Yixiang Gan and Dr. Dorian Hanaor. 
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Stress-Dependent Electrical Contact 
Resistance at Fractal Rough Surfaces 


Chongpu Zhai'; Dorian Hanaor*; Gwénaélle Proust®; and Yixiang Gan* 


Abstract: The electrical contact resistance between contacting rough surfaces was studied under various compressive stresses. The samples 
considered here were isotropically roughened aluminium disks with upper and lower surfaces modified through polishing and sand blasting 
using different sized glass beads. Fractal geometry and roughness descriptors, including root mean square values of roughness and slope, were 
used to describe the topography of sample surfaces, based on the digitized profiles obtained from interferometry-based profilometry. The 
electrical contact resistances at the interfaces were obtained by applying a controlled current and measuring the resulting voltage, through the 
following scenarios: (1) over time for various applied testing currents, the resistance relaxation curves were measured at constant loads; 
(2) through voltage-current characteristics by means of a logarithmic sweeping current, the influence of the testing current on the electrical 
response of contacting rough surfaces was evaluated; and (3) for a given testing current, the electrical resistance through interfaces of different 
surface structures was measured under increasing compressive stresses. The experimental results show that the measured resistance depends 
closely on the measurement time, testing current, surface topology, and mechanical loading. At stresses from 0.03 to 1.18 MPa, the electrical 
resistance as a function of applied normal stress is found to follow a power law relation, the exponent of which is closely linked to the surface 


topology. DOI: 10.1061/(ASCE)EM.1943-7889.0000967. © 2015 American Society of Civil Engineers. 


Author keywords: Rough surfaces; Electrical contact resistance; Branly effect; Fractal dimension. 


Introduction 


The properties of the electrical contact resistance (ECR) at electri- 
cal connections are of tremendous importance in many engineering 
applications, including resistance spot welding, diagnostic tribol- 
ogy, signal, and current transmission (Crinon and Evans 1998; 
Kogut and Komvopoulos 2003; Kogut and Komvopoulos 2005; 
Slade 2013). For most applications, e.g., electronic connectors 
(Bryant and Jin 1991), switches, and electrode structures of bat- 
teries (Meulenberg et al. 2003), a low and stable contact resistance 
is sought after. In these situations, numerous factors can affect the 
current flow through contacting surfaces, e.g., surface topography, 
environmental conditions, mechanical loads, the presence of a coat- 
ing layer, and applied voltage. Vibrations, fretting, thermal shock, 
and chemical contamination are principal failure mechanisms at 
electrical contacts often resulting in the malfunction of electric cir- 
cuits. Effective electrical contacts require reliable mechanical con- 
tact in order to avoid failure. Surface characteristics, such as 
roughness, curvature, and fractal dimension, play an important role 
in facilitating the closing and opening of electric circuits created by 
the close contact between surfaces (Bryant and Jin 1991; Oh et al. 
1999; Kogut and Etsion 2000; Falcon et al. 2004; Kogut 2005). 
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The true area of contact at an interface is considerably smaller 
than the apparent or nominal contact area because of the existence 
of surface roughness (Archard 1957; Bowden and Williamson 
1958; Greenwood and Williamson 1958). When two rough surfa- 
ces are squeezed together contact is made through individual asper- 
ities with contact patches that can extend in size down to the 
nanoscale. These contact junctions exhibit electrical and mechani- 
cal properties that may diverge from bulk properties (Jackson et al. 
2012). Current flowing through rough interfaces is scattered across 
many contacting asperities with electronic transport involving 
multiple mechanisms, including quantum tunneling (Yanson 
et al. 1998; Foley et al. 1999; Agrait et al. 2003), Sharvin contact 
(Sharvin 1965), and Holm contact (Holm and Holm 1967), depend- 
ing on the size of contacting junctions and the mean free path of 
electrons. Early work by Holm and Holm (1967) concluded that the 
ECR is affected by both the constriction resistance resulting from 
the limited areas of true contact at an interface and interfacial 
resistance due to the inevitable presence of resistive surface films, 
such as oxide layers. When two metal surfaces are compressed 
together with sufficient pressure, surface asperities can penetrate 
the oxide layer thus forming metal-to-metal contact patches. When 
the size of contacting asperities becomes larger than the mean free 
path of electrons, Holm contact will be the dominant transport 
mechanism, resulting in a relative low resistance. This concept 
has been further developed by Chang et al. (1987), Ciavarella et al. 
(2004, 2008), and Jackson et al. (2009). These theoretical and com- 
putational studies assumed that individual microcontacts formed 
by asperities, where the circuit continuity was established, gov- 
erned specific electron transport regions, limited by the measuring 
resolution. 

The asperities of contacting surfaces tend to exhibit complex 
geometries and structures at a wide range of length scales, 
governing physical properties and interfacial phenomena and giv- 
ing rise to constriction resistance. The fractal topography based 
on scale-invariant parameters provides an effective means for 
modeling engineering surfaces with random self-affine multiscale 
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properties (Yan and Komvopoulos 1998; Kogut and Komvopoulos 
2004). A general ECR theory based on fractal geometry was pro- 
posed in order to describe the effects of contact loads, elastic-plastic 
deformation of the contacting asperities, surface topography, and 
material properties on size-dependent electrical resistance of the 
microcontacts comprising the true contact area (Mikrajuddin et al. 
1999; Kogut and Komvopoulos 2003). Kogut et al. (2005) further 
investigated conductive rough surfaces separated by a thin insulat- 
ing film. Semiempirical power-law type correlations between the 
contact resistance and the normal pressure have also been pro- 
posed. For contacting surfaces separated by superficial oxide or 
impurity layers, the power law is found, with the exponent value 
possibly greater than 1 (Milanez et al. 2003; Falcon and Castaing 
2005; Paggi and Barber 2011). The power-law relation between the 
applied normal load and the contact conductance at rough surfaces 
has been reported on the basis of the incremental stiffness, which 
was found to be linearly proportional to the contact conductance 
(Barber 2003; Pohrt and Popov 2012; Pohrt and Popov 2013; 
Hanoar et al. 2015). 

Furthermore, the influencing factors of ECR include the physi- 
cal and chemical origins and details of the oxidation processes in 
corrosive environment (Sun et al. 1999; Sun 2001), surface diffu- 
sion (Crinon and Evans 1998; Ogumi and Inaba 1998) and fretting 
corrosion (Bryant 1994). Thin oxide or hydroxide layers, acting as 
resistive films with typical thickness ranging from | to 10 nm, tend 
to form at metallic surfaces, covering the contacting surface, and 
increasing the contact resistance (Bryant and Jin 1991). In conduc- 
tion through interfaces, insulating layers can bring about high 
electrical resistance under conditions of low current. With an 
increasing current this high initial resistance decreases by several 
orders of magnitude showing nonlinear conductivity phenomena, 
which is called the Coherer effect or Branly effect. The process 
is featured by voltage creep, hysteresis loops, and voltage satura- 
tion effects (Castaing and Laroche 2004; Falcon and Castaing 
2005; Bourbatache et al. 2012; Tekaya et al. 2012). Many mech- 
anisms of resistance modification have been reported 1n conduction 
phenomena through rough interfaces, including electron tunneling 
through oxide layers and voids (Creyssels et al. 2007), dielectric 
breakdown of oxide layers (Dorbolo et al. 2002), localized current- 
induced welding (Falcon and Castaing 1993), percolation collective 
process (Falcon and Castaing 2005), and chemical disorder arising 
with random composition (Creyssels et al. 2009). 

Despite recent progress of nanoscale testing technology and 
surface morphology characterization, a stress-dependence of elec- 
trical conduction behavior through rough interfaces with random 
multiscale texture remains largely unknown. Particularly, existing 


experimental results are limited. In this paper, the measured ECR 
of aluminium disk stacks loaded in compression is presented. The 
surfaces of the disks were modified by polishing and sand blasting 
in order to obtain a range of rough surfaces. Compressive loads, 
testing current, and time were varied to study their influence on 
the ECR. 


Sample Preparation and Surface Characterization 


In this paper, round disks, with a diameter of 25 mm, made of 
aluminium alloy were used as specimens. Surfaces were polished 
followed by sand abrasive blasting processes using different sized 
particles to modify the surface details and structures at various 
scales (Hanaor et al. 2013). For each individual sample, both 
top and bottom were equally treated through standard polishing 
and sand blasting procedures. The average diameters of the two 
selected groups of glass beads used in blasting treatments were 
50 and 300 ym. Fig. 1 shows scanning electron microscope (SEM) 
images of the different surface types used in this work. Samples that 
have been blasted using glass beads of 50 zm are found to exhibit 
the most complex and roughest surface texture. The treated alumin- 
ium surfaces were scanned using an optical surface profilometer 
(NanoMap 1OOOWLI, AEP Technology, California) to obtain 
three-dimensional digitized topographies. Subsequently, RMS 
roughness, RMS slopes, and fractal dimension were utilized to 
characterize and compare the surface geometries. 

As shown in Table 1, three types of samples with distinct surface 
details were prepared and characterized prior to electrical testing, 
with Sl, S2, and S3 representing the polished samples, samples 
blasted using 300 zm diameter glass beads, and samples blasted 
using 50 zm diameter glass beads, respectively. Before the standard 
sand blasting treatment, all the sample surfaces were polished and 
prepared using several polishing steps with the final step using 
1 um diamond suspension. Then, the processed samples were prop- 
erly cleaned by water and compressed air to remove any embedded 
glass beads. Cleaning with ethanol and heat treatment (from 110 to 
120°C) were also applied to remove surface contamination and 
moisture. For each surface, the mean values of roughness descrip- 
tors with standard deviations over ten | x 1 mm scans from differ- 
ent samples are shown in Table |. The scaled triangulation method 
(Dubuc et al. 1989; De Santis et al. 1997; Zahn and Zosch 1999) 
was used for the calculation of fractal dimension values. A power 
law relationship is found between the calculated surface area and 
the length resolution of digitized scan of all the surfaces, exhibiting 
self-affinity over a range of length scales (from | to 100 pm). 





Fig. 1. SEM images of aluminium samples with different surface treatments: (a) polishing treatment; (b) blasted with 300 jzm diameter glass beads; 


(c) blasted with 50 ym diameter glass beads 
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Table 1. Sample Surface Characterization with Different Treatments 


Sample RMS roughness RMS slope Fractal 

type Surface treatment Rems/um Rs dimension D/triangulation 
S1 Polished 0.057 + 0.005 0.009 + 0.001 2.093 + 0.0620 

S2 Blasted by glass beads of 300 jum diameter 4.179 + 0.194 0.224 + 0.015 2.051 40.0217 

S3 Blasted by glass beads of 50 wm diameter 2.970 + 0.276 0,202 = :0,010 2.626 + 0.0174 


Surfaces blasted with glass beads of 300 um diameter reveal 
the highest RMS slope and RMS roughness. The lowest values 
of RMS roughness, RMS slope and fractal dimension are found 
for the polished surfaces. The values of lower-roll-off wavelength 
in the power spectrum of all of the surfaces used were found to 
be ~100 ym, which is remarkably smaller than the disk diameter 
(25 mm). 


Results and Discussion 


Experimental Setup 


For each type of sample with similar surface characteristics, elec- 
trical resistances were measured for stacks of 11 samples, giving 10 
rough to rough interfaces, by means of a source/measurement unit 
(Agilent B2902A, Keysight Technologies), under various compres- 
sive loading forces, as shown in Fig. 2. In this experimental setting, 
the resistance formed by 10 specimen interfaces placed between 
two polished plates of the same material instead of one single inter- 
face is measured, aiming to achieve a higher precision, larger linear 
range, and better robustness against measurement noises from the 
connecting wires, loading device, and measurement unit. A stack of 
samples can also diminish effectively the experimental errors from 
the variation of samples with the same surface treatment. Before the 
experiment, the sample stack was aligned by rotating the top pres- 
sure head of the loading machine and using a piece of rubber placed 
between the pressure head and the top polished plate. 

Prior to the measurement of stress-dependent resistance, the 
following tests were performed: (1) resistance creep tests and 
(2) sweeping current tests, to exclude the influences from the ap- 
plied current and measurement time. In the following sections, 
the results of these tests will be presented and discussed. Finally, 


Current through 
rough interface 





—> Holm contact 


a> Sharvin effect 
—>p> Tunneling through \ 
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ayy Tunneling through 
the void 


the stress-dependent ECR was measured using the applied current 
of 10 mA and measuring time of 0.01 s at each individual stress 
level. 


Resistance Creep Tests 


The relaxation dynamics of resistance for a given constant current 
were first investigated to find out the effect of current on the 
measured resistance over time. The resistance degradation curves 
at various applied testing currents of 10, 20, 80, 320, and 1,280 mA 
were recorded, under an applied pressure of 0.061 MPa, as shown 
in Fig. 3(a). Then, a constant current (320 mA) was applied to 
study the resistance relaxation under various applied compressive 
stresses, ranging from 0.031 to 0.490 MPa, as shown in Fig. 3(b). 
Multiple tests were performed for each individual loading condi- 
tion, but for clarity here not all the experimental data collected were 
shown and the trends shown in Fig. 3 are true for all the data col- 
lected. The measured resistance under various loads and current 
conditions gradually decreases with respect to time. For tests 
carried out under a constant normal load, shown in Fig. 3(a), a 
higher current brings about a more significant decrease in measured 
resistance. More specifically, the highest current of 1,280 mA over 
200 s was found to cause a decrease of around 15% with reference 
to the initial measured resistance, whereas a variation of less than 
1% was achieved at the lowest current (10 mA). For the tests with 
the same testing current (320 mA), a higher pressure reduces the 
trend, which is presented in Fig. 3(b). The measured resistance of 
the samples under an applied pressure of 0.490 MPa exhibited a 
slight decrease of less than 2% with reference to the initial resis- 
tance, whereas a drop of approximately 15% was found under an 
applied pressure of 0.031 MPa. For all the tests when the conduc- 
tion time is less than 1 s, the decline of the measured resistance due 
to the current is less than 0.1%. 


Load 
Bubber Resistance 
Measurement 






a 


Fig. 2. Experimental setup for the resistance measurement of a stack of rough surfaces and schematic of the current flow through rough interfaces by 


means of Holm contact, Sharvin contact and quantum tunneling 
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Fig. 3. Typical resistance relaxation curves over 200 s, with R; being the initial measured resistance: (a) with constant applied stress being 0.061 MPa 
under various testing current (10, 20, 80, 320, and 1,280 mA with the initial resistance being 1.57, 1.46, 0.93, 0.76, and 0.23 (2, respectively); (b) with 
the constant testing current being 320 mA under various compressions (0.031, 0.061, 0.122, and 0.490 MPa with the initial resistance of 2.01, 0.76, 


0.59, and 0.18 (2, respectively) 


Sweeping Current Tests 


The influence of the testing current on measured ECR was assessed 
by sweeping current tests. Here the applied sweep current test con- 
sisted of two phases with the /oading phase (P1) having a logarith- 
mically increasing current from 0.0001 to 1.5 A and the unloading 
phase (P2) having a logarithmically decreasing current from 1.5 A 
back to 0.0001 A. During the sweep procedure, the voltage was 
recorded at a data acquisition frequency of 2 kHz corresponding 
to the electrical current imposed. The normal pressure applied 
on the measured samples was 30 N, corresponding to a stress level 
of 0.061 MPa. The sweep process, including both phases, was ac- 
complished within 0.2 s in order to avoid the occurrence of elec- 
trical degradation over time. As discussed in the previous section, 
when the conduction time is shorter than | s, the time-dependant 
reduction in resistance is negligible. Figs. 4(a and b) show the mea- 
sured voltage and contact resistance with respect to testing current 
for the three types of samples. The shape of voltage hysteresis loops 
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is seen to be consistent for the three surface types while presenting 
different voltage levels. In Fig. 4(b), the measured results of all 
three surfaces demonstrate similar trends known as the Branly ef- 
fect (Castaing and Laroche 2004; Falcon and Castaing 2005), 
1.e., the measured resistance begins to drop after the testing current 
reaches a certain value, after which the resistance appear to be irre- 
versible as the current sweeps back. The corresponding threshold 
current values for S1, S2, and S3 are around 150, 80, and 60 mA, 
respectively, and the values seem to inversely correlate with the 
fractal dimension of the surface. However, the definite correlation 
between the RMS roughness and the threshold current can only be 
concluded by extending the types of surfaces. The use of RMS 
roughness is limited in predicting interface behaviors, whereas 
the fractal dimension, a cross-scale descriptor, is extensively used 
to interpret surface phenomena, such as conduction properties 
at rough interfaces (Kogut and Komvopoulos 2005; Persson 
2006). The transition at the threshold current may come from the 
electrothermal coupling of the microcontacts at rough interfaces. 
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Fig. 4. Typical measured results for three types of samples using sweeping currents under a constant stress of 0.061 MPa; the data presented by the 
solid lines describes the first phases with increasing testing currents and the dashed lines show the second phases with decreasing currents: (a) variation 


of the measured voltage; (b) resistance-current characteristics 
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For testing currents higher than approximately 5 mA and lower 
than the threshold current values, the measured resistances remain 
stable on two plateaus in both P| and P2, with a larger value ob- 
tained for P1, shown in Fig. 4(b). In other words, ohmic behavior 
at low levels of electrical current was found for all three types of 
surfaces. As shown in Fig. 4(b), at low testing currents (lower than 
1 mA) in either Pl! or P2, the measured resistances exhibit insta- 
bility, especially for polished samples. Measurement noise for all 
three surfaces is observed at similar levels, which are amplified in 
logarithmic scale, being prominent when the measured resistances 
are small. 

For each type of surfaces, the range of measured voltage and 
resistance was found to vary in five different tests under the same 
experimental conditions, e.g., the measured resistance of polished 
samples spreads from 0.5 to 3 (2. From the experimental results, the 
surfaces described by higher fractal dimension demonstrate larger 
values of measured resistances. Specifically, the smallest resistan- 
ces are observed from polished samples and the largest resistances 
are measured from the samples blasted using glass beads of 50 wm 
diameter. The polished samples present the smoothest and the most 
stable trends among the three types of surfaces. High level noises 
are found for surfaces after sand blasting procedures showing more 
complex topographies, which may cause the sensitivity to the mi- 
crovibration from the loading device and the electrical noise from 
the measurement unit. 

By rearranging the sample stacking order or by rotating the sam- 
ples to achieve different contact configurations one can retrieve the 
voltage hysteresis loops and repeat the processes depicted in Fig. 4. 
Despite this, all the samples used were repolished and sand blasted 
again between successive tests to avoid the accumulation of surface 
modification by current flow as the sweep procedure may bring 
about some localized modification of the surface characteristics 
or possible changes of surface chemical composition. 

A further test was done by driving successive back-and-forth 
scanning current cycles with increasing current range. A typical 
obtained result for polished samples is shown in Fig. 5. Here each 
successive current cycle was of 0.2 s duration and compressive 
pressure was maintained at 0.02 MPa during all nine sweep cycles. 
The sweeping current range was increased by a factor of /2 at each 
consecutive cycle until the current sweep reached 1.5 A. There 
exists a pause of around 10 s for data exporting and setting of 


Measured voltage (V) 
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(a) Testing current (A) 


the source/measurement unit. Additionally, tests with a longer 
pause up to several minutes were conducted and presented the same 
trends. For all nine cycles, the electrical current flows at different 
levels can modify the ECR of interfaces at different degrees. At 
each individual cycle, the measured voltages and resistances tend 
to be reversible when the testing current is lower than the maximum 
current of the previous cycle, 1.e., the measured voltages and resis- 
tances in PI branch are inversely along the track of P2 of the pre- 
vious cycle, exhibiting ohmic conduction properties. The envelope 
line constituted by P1 of the first cycle, P2 of the last cycle, and the 
resistance declining segments in P1 of all nine cycles, exhibits a 
similar curve shapes as that of each individual cycle. The threshold 
current values at the turning point for measured resistance in all 
cycles grow from around 100 mA to 1 A, along with the increasing 
range of the imposed sweeping current. 

Even though all the tests were performed over a short duration, 
1.e., 0.2 s, the Joule heating generated from the testing current may 
still have an effect on the measured resistance. In Fig. 5(b), evident 
upward and downward trends are observed at the end of P1 and at 
the beginning of P2, respectively. The aluminium has a lower heat 
capacity than the oxide. The Joule heating is likely to contribute to 
the trends to a certain extent. The irreversible Joule heating from 
the injected current, which accumulates with respect to time, cannot 
fully explain the reversible upward and downward trends at the end 
of Pl and the beginning of P2. Another possible reason can be the 
charging and discharging processes occurring during the tests. The 
contact of rough interfaces can be regarded as a complex network 
of resistors and capacitors that vary with the applied pressure and 
injected current. Moreover, this process is very similar to the con- 
duction behaviors in semiconducting devices showing nonlinear 
and reversible properties. 

A further experimental study was conducted in order to 
shed light on relations between the applied load and ECR. For pol- 
ished surfaces, five typical voltage-current loops and resistance- 
current characteristics under various loads ranging from 0.031 to 
0.490 MPa were obtained and are shown in Fig. 6. Under a high 
compressive pressure (more than 0.490 MPa), the loops become 
flat, i.e., P2 follows the same path as Pl. The Branly effect tends 
to be harder to capture when the pressed surfaces were in contact 
at sufficiently high stress values. A high level of applied stress 
leads to better stability and repeatability of ECR measurements. 


NO 





Measured resistance (Q) 
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Fig. 5. Typical measured results for a stack of polished samples with nine successive back-and-forth sweeping current cycles with range increasing by 
a factor of \/2 from 0.1024 to 1.5000 A (0.1024, 0.1448, 0.2048, 0.2896, 0.4096, 0.5793, 0.8192, 1.1585 and 1.5000 A) under a constant stress of 
0.02 MPa: (a) hysteresis loops of voltage; (b) resistance-current characteristics 
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Fig. 6. Typical measured results for polished samples using current sweep under various stresses with solid lines expressing the first phase and the 
dashed line showing the second phase: (a) variation of the measured voltage; (b) resistance-current characteristics 


The adjacent asperities may move closer to each other under a 
higher contact pressure, which results in the rise of the contacting 
area, thus enhances the conduction. The enlarged contacting area 
can carry more current without the microsoldering of the contact 
points. In engineering practice, the contact resistance is usually re- 
duced by squeezing the contacting components together with 
greater force, which can bring about a larger true contact area. With 
the sweep tests shown here, the electrical current seems to be also 
capable of broadening the current path. 


Stress-Dependent Electrical Contact Resistance 


The results shown in Figs. 3—6 indicate a likely change of the 
contact network between the surfaces, where the contacting asper- 
ities can be regarded as a network of resistors and capacitors chang- 
ing with applied current, mechanical load, and time. The true 
contact area increases linearly with the applied compression, result- 
ing in improved conduction (Kogut 2005; Kogut and Komvopoulos 
2005). Meanwhile, the electrical current results in the physical and 
chemical modification of sample surfaces, which involves many 
processes, including the rupture of the oxide layer due to compres- 
sion and the localized heating induced by current. 

With the tests presented above, the conclusions are (1) for the 
studied systems, when the imposed current 1s less than 50 mA, im- 
pact from the current on the measured resistance is negligible, even 
at low levels of applied stress (lower than 0.020 MPa). However, a 
small testing current of less than | mA can lead to a high-level of 
measurement noise; and (2) when the conduction time of electrical 
current (less than 50 mA) is less than 1 s, the influence of the cur- 
rent on resistance can be ignored, especially in cases of high pres- 
sures (higher than 0.5 MPa). 

For samples exhibiting different surface morphologies, the elec- 
trical resistances were measured under various stresses, as shown in 
Fig. 7. For each type of sample, five series of tests were conducted 
and the resistances were evaluated at 16 different stress levels from 
0.020 to 8.936 MPa. The testing current was set as 10 mA and the 
measured time was 0.01 s for each individual data point in order to 
minimize the influence of the testing current on the measurement. 

As shown in Fig. 7, the measured resistances of disk stacks with 
different surface features decrease considerably as the compressive 
pressure is increased, converging to a value close to the bulk resis- 
tance of the material. For the identical loading stress level, samples 
blasted with 50 yum sized glass beads usually present the highest 
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resistance among all three types of samples. At low levels of 
applied stress (less than 0.5 MPa), the measured resistance 1s spread 
across a wider range. A few groups of the testing results were not 
included when jumps occur in the resistance measurement under 
constant or smoothly increasing loads. The unexpected vibration, 
electrical noise and some limitations of the experimental setting can 
contribute to the measurement uncertainty. By fitting the resistance/ 
pressure curves from 0.031 to 1.176 MPa it is observed that 
the measured contact resistance is a power function of the compres- 
sive stress at certain range of loading with the exponents being 
—0.816, —1.026, and —1.494, respectively, for polished surface 
(S1), blasted using 300 um (S2), and 50 ym (S3) diameter glass 
beads. The absolute values of power exponents demonstrate an in- 
creasing trend with the rise of fractal dimension indicated in Table 1. 
When a loading force of 4,384 N is applied, corresponding to a 
stress level of 8.936 MPa, the measured resistance can be as small 
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Fig. 7. Stress-dependent ECR of different surfaces under various 
loading levels with a testing current of 10 mA, with EF being the value 
of the Young’s modulus of the tested material, A being the projected 
area of the tested samples, Rp being 0.06 22, including the combined 
resistance of bulk material of identical size as the disk stack (equaling 
2.53 Q), wires and connections used in the experimental setting, and 
R, being 1 Q 
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as 0.0578, 0.0271, and 0.0303 (2, corresponding to S1, S2, and S3. 
The resistance measured at high loads for the three types of tested 
surfaces appears to be inversely correlated with the values of RMS 
roughness and RMS slope. Further experiments with more types of 
surface are necessary to reach a general conclusion on the relation- 
ship between the stress-dependent electrical resistance and RMS 
roughness and RMS slope. 


Conclusions 


An experimental investigation was performed on the electrical con- 
tact resistance of rough surfaces. The experimental results show 
that the measured contact resistances of a stack of rough aluminium 
samples rely strongly on surface topology, mechanical loading, as 
well as testing current and time. The typical resistance creep curves 
over 200 s were recorded to investigate the current effect on the 
measured resistance over time. Tests with sweeping current under 
various compressive loads were carried out to study the influence of 
the current on the electrical contact resistance. For low-level elec- 
trical currents of the orders of a few mA, the electrical resistance is 
constant, exhibiting a linear ohmic behavior whereas the rise of the 
electrical current brings about a decrease of the electrical resistance. 
An increased compressive load results in the weakening of the ef- 
fects of testing current on the measured resistance and decreases the 
time dependent resistance relaxation. Similarly, a higher testing 
current can also cause the reduction of the impact from the loading 
force. With testing current set as 10 mA and the measuring time 
being 0.01 s for each measured resistance, the impact of the current 
can be ignored for all three types of samples under varying loading 
levels. Under this condition, the measured resistances of stacks of 
samples decrease continuously under increasing stresses, approach- 
ing the resistance for bulk material at high loads. The results also 
demonstrate that a power law relationship exists between the mea- 
sured electrical resistance and normal stress across a certain stress 
range, with the absolute values of exponents rising with the fractal 
dimension of the surfaces. 


Acknowledgments 


Financial support for this research from the Australian Research 
Council through grants DE130101639 and Civil Engineering 
Research Development Scheme (CERDS) in the School of Civil 
Engineering at The University of Sydney is greatly appreciated. 


References 


Agrait, N., Yeyati, A. L., and Van Ruitenbeek, J. M. (2003). “Quantum 
properties of atomic-sized conductors.” Phys. Rep., 377(2), 81-279. 

Archard, J. (1957). “Elastic deformation and the laws of friction.” Proc. R. 
Soc. London, Ser. A, 243(1233), 190-205. 

Barber, J. (2003). “Bounds on the electrical resistance between contacting 
elastic rough bodies.” Proc. R. Soc. London, Ser. A, 459(2029), 53-66. 

Bourbatache, K., Guessasma, M., Bellenger, E., Bourny, V., and Tekaya, A. 
(2012). “Discrete modelling of electrical transfer in multi-contact sys- 
tems.” Granular Matter, 14(1), 1-10. 

Bowden, F., and Williamson, J. (1958). “Electrical conduction in solids. I: 
Influence of the passage of current on the contact between solids.” Proc. 
R. Soc. London, Ser. A, 246(1244), 1-12. 

Bryant, M. D. (1994). “Resistance buildup in electrical connectors due to 
fretting corrosion of rough surfaces.” Compon. Packag. Manuf. Tech- 
nol. Part A IEEE Trans., 17(1), 86—95. 

Bryant, M. D., and Jin, M. (1991). “Time-wise increases in contact resis- 
tance due to surface roughness and corrosion.’ Compon. Hybrids 
Manuf. Technol. IEEE Trans., 14(1), 79-89. 


© ASCE 


B4015001-7 


Castaing, B., and Laroche, C. (2004). “Turbulent” electrical transport in 
copper powders.” Europhys. Lett., 65(2), 186-192. 

Chang, W., Etsion, I., and Bogy, D. B. (1987). “An elastic-plastic model for 
the contact of rough surfaces.” J. Tribol., 109(2), 257-263. 

Ciavarella, M., Dibello, S., and Demelio, G. (2008). “Conductance of 
rough random profiles.” Int. J. Solids Struct., 45(3), 879-893. 

Ciavarella, M., Murolo, G., and Demelio, G. (2004). “The electrical/ 
thermal conductance of rough surfaces—the Weierstrass-Archard multi- 
scale model.” Int. J. Solids Struct., 41(15), 4107-4120. 

Creyssels, M., Dorbolo, S., Merlen, A., Laroche, C., Castaing, B., and 
Falcon, E. (2007). “Some aspects of electrical conduction in granular 
systems of various dimensions.” Eur. Phys. J. E Soft Matter Biol. Phys., 
23(3), 255-264. 

Creyssels, M., Falcon, E., and Castaing, B. (2009). “Experiment and theory 
of the electrical conductivity of a compressed granular metal.” Proc., 
6th Int. Conf. on the Micromechanics of Granular Media, American 
Institute of Physics, Melville, NY, 123-126. 

Crinon, E., and Evans, J. (1998). “The effect of surface roughness, oxide 
film thickness and interfacial sliding on the electrical contact resistance 
of aluminium.” Mater. Sci. Eng. A, 242(1), 121-128. 

De Santis, A., Fedi, M., and Quarta, T. (1997). “A revisitation of the tri- 
angular prism surface area method for estimating the fractal dimension 
of fractal surfaces.” Ann. Geophys., 40(4), 811-821. 

Dorbolo, S., Ausloos, M., and Vandewalle, N. (2002). “Hysteretic behavior 
in metallic granular matter.’ Appl. Phys. Lett., 81(5), 936-938. 

Dubuc, B., Zucker, S., Tricot, C., Quiniou, J., and Wehbi, D. (1989). 
“Evaluating the fractal dimension of surfaces.” Proc. R. Soc. London 
A. Math. Phys. Sci., 425(1868), 113-127. 

Falcon, E., and Castaing, B. (1993). “Electrical properties of granular mat- 
ter: From” Branly effect” to intermittency.” Proc., Geotechnical Man- 
agement of Waste and Contamination, CRC Press, London, U.K., 323. 

Falcon, E., and Castaing, B. (2005). “Electrical conductivity in granular 
media and Branly’s coherer: A simple experiment.” Am. J. Phys., 
73(4), 302-307. 

Falcon, E., Castaing, B., and Creyssels, M. (2004). “Nonlinear electrical 
conductivity in a ID granular medium.” Eur. Phys. J. B, 38(3), 
475-483. 

Foley, E., Candela, D., Martini, K., and Tuominen, M. (1999). “An under- 
graduate laboratory experiment on quantized conductance in nanocon- 
tacts.” Am. J. Phys., 67(5), 389-393. 

Greenwood, J., and Williamson, J. (1958). “Electrical conduction in solids. 
Il. Theory of temperature-dependent conductors.’ Proc. R. Soc. 
London, Ser. A, 246(1244), 13-31. 

Hanaor, D., Gan, Y., and Einav, I. (2013). “Effects of surface structure 
deformation on static friction at fractal interfaces.” Géotechnique Lett., 
3(2), 52-58. 

Hanaor, D., Gan, Y., and Einav, I. (2015). “Contact mechanics of fractal 
surfaces by spline assisted discretisation.” Int. J. Solids Struct., 
59, 121-131. 

Holm, R., and Holm, E. A. (1967). Electric contacts, Springer, New York. 

Jackson, R. L., Crandall, E. R., and Bozack, M. (2012). “An analysis of 
scale dependent and quantum effects on electrical contact resistance be- 
tween rough surfaces.” Proc., Electrical Contacts (Holm), 2012 [EEE 
S&th Holm Conf. on, IEEE, New York, 1-10. 

Jackson, R. L., Malucci, R. D., Angadi, S., and Polchow, J. R. (2009). “A 
simplified model of multiscale electrical contact resistance and com- 
parison to existing closed form models.” Proc., Electrical Contacts, 
2009 Proc., 55th IEEE Holm Conf. on, IEEE, New York, 28-35. 

Kogut, L. (2005). “Electrical performance of contaminated rough surfaces 
in contact.” J. Appl. Phys., 97(10), 103723. 

Kogut, L., and Etsion, [. (2000). “Electrical conductivity and friction 
force estimation in compliant electrical connectors.” Tribol. Trans., 
43(4), 816-822. 

Kogut, L., and Komvopoulos, K. (2003). “Electrical contact resistance 
theory for conductive rough surfaces.” J. Appl. Phys., 94(5), 
3153-3162. 

Kogut, L., and Komvopoulos, K. (2004). “Electrical contact resistance 
theory for conductive rough surfaces separated by a thin insulating 
film.” J. Appl. Phys., 95(2), 576-585. 


J. Eng. Mech. 


J. Eng. Mech., B4015001 


Downloaded from ascelibrary.org by UNIVERSITY OF SYDNEY on 10/01/16. Copyright ASCE. For personal use only; all rights reserved. 


Kogut, L., and Komvopoulos, K. (2005). “Analytical current-voltage rela- 
tionships for electron tunneling across rough interfaces.” J. Appl. Phys., 
Oitt), OF 2701, 

Meulenberg, W., Uhlenbruck, S., Wessel, E., Buchkremer, H., and St6ver, 
D. (2003). “Oxidation behaviour of ferrous alloys used as interconnect- 
ing material in solid oxide fuel cells.” J. Mater. Sci., 38(3), 507-513. 

Mikrajuddin, A., Shi, F. G., Kim, H., and Okuyama, K. (1999). “Size- 
dependent electrical constriction resistance for contacts of arbitrary size: 
From Sharvin to Holm limits.” Mater. Sci. Semicond. Process., 2(4), 
321-327. 

Milanez, F. H., Yovanovich, M. M., and Culham, J. R. (2003). “Effect of 
surface asperity truncation on thermal contact conductance.” Compon. 
Packag. Technol. IEEE Trans., 26(1), 48—54. 

Ogumi, Z., and Inaba, M. (1998). “Electrochemical Lithium Intercalation 
within carbonaceous materials: Intercalation processes, surface film 
formation, and lithium diffusion.” Bull. Chem. Soc. Japan, 71(3), 
521-534. 

Oh, H.-J., Jang, K.-W., and Chi, C.-S. (1999). “Impedance characteristics 
of oxide layers on aluminium.” Bull. Korean Chem. Soc., 20(11), 
1340-1344. 

Paggi, M., and Barber, J. (2011). “Contact conductance of rough surfaces 
composed of modified RMD patches.” Int. J. Heat Mass Transf., 
54(21), 46644672. 

Persson, B. N. J. (2006). “Contact mechanics for randomly rough surfaces.” 
Surf. Sci. Rep., 61(4), 201-227. 


© ASCE 


B4015001-8 


Pohrt, R., and Popov, V. L. (2012). “Normal contact stiffness of elastic 
solids with fractal rough surfaces.” Phys. Rev. Lett., 108(10), 104301. 

Pohrt, R., and Popov, V. L. (2013). “Contact stiffness of randomly rough 
surfaces.” Sci. Rep., 3, 3293. 

Sharvin, Y. V. (1965). “On the possible method for studying Fermi surfa- 
ces.” Zh. Eksperim. i Teor. Fiz., 48, 984-985. 

Slade, P. G. (2013). Electrical contacts: Principles and applications, CRC 
Press, London, U.K. 

Sun, M. (2001). “Conductivity of conductive polymer for flip chip bonding 
and BGA socket.” Microelectron. J., 32(3), 197-203. 

Sun, M., Pecht, M. G., Natishan, M. A. E., and Martens, R. I. (1999). “Life- 
time resistance model of bare metal electrical contacts.” Adv. Packag. 
IEEE Trans., 22(1), 60-67. 

Tekaya, A., Bouzerar, R., and Bourny, V. (2012). “Influence of surface 
topology on the electrical response of many bead assemblies.” A/P 
Adv., 2(3), 032108. 

Yan, W., and Komvopoulos, K. (1998). “Contact analysis of elastic-plastic 
fractal surfaces.” J. Appl. Phys., 84(7), 3617-3624. 

Yanson, A., Bollinger, G. R., Van den Brom, H., Agrait, N., and 
Van Ruitenbeek, J. (1998). “Formation and manipulation of a metallic 
wire of single gold atoms.” Nature, 395(6704), 783-785. 

Zahn, W., and Zodsch, A. (1999). “The dependence of fractal dimension 
on measuring conditions of scanning probe microscopy.” Fresenius’ 
J. Anal. Chem., 365(1-3), 168-172. 


J. Eng. Mech. 


J. Eng. Mech., B4015001 


4.2 Paper 4: Interfacial electro-mechanical behaviour 


The inter-particle electrical property is playing a significant role in governing the electrical transport 
of granular materials. This paper presented here studies interfacial stiffness and electrical contact 
resistance, which are both closely related to real contact area, allow us to understand how two 
contacting solids with rough surfaces interact with each other over a wide range of length scales. The 
results of our original experiments show that the interfacial contact stiffness and electrical 
conductance follow power-law functions with respect to the normal force over a wide range of applied 
loads. However, different values of power-law exponents were found. Specifically, the power-law 
exponent for electrical contact resistance is approximately double of that for normal interfacial 
stiffness. This is because the transport behaviour at the rough contact at macro scales is dominated by 
the contact behaviour at lower scales. Moreover, the microstructures represented by surface roughness 
bring about strong interconnections of electrical-mechanical coupling. This work provides a first- 
order investigation connecting interfacial mechanical and electrical behaviour, applicable to studies of 
electrical transport in granular materials with numerous contacts. Both the interfacial stiffness 
presented in Chapter 3 and electrical contact resistance shown in this chapter are fundamentally 
determined by the real contact area at rough interface. By relating the two parameters, one can obtain 
a comprehensive understanding of the rough contact and the electrical transport at rough interfaces. 
The information obtained in Chapter 3 and 4 will be scaled upwards to granular systems at macroscale 
through contact network built by individual inter-particle contacts to interpret the stress-dependent 
electrical conduction in granular materials, shown in the following two chapters. 


This paper has been published in Extreme Mechanics Letters in 2016. I was the primary researcher 
and author of this paper, being supervised by Dr. Yixiang Gan and Dr. Dorian Hanaor. 
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1. Introduction 


ABSTRACT 


In a range of energy systems, interfacial characteristics at the finest length scales 
strongly impact overall system performance, including cycle life, electrical power loss, 
and storage capacity. In this letter, we experimentally investigate the influence of surface 
topology on interfacial electro-mechanical properties, including contact stiffness and 
electrical conductance at rough surfaces under varying compressive stresses. We consider 
different rough surfaces modified through polishing and/or sand blasting. The measured 
normal contact stiffness, obtained through nanoindentation employing a partial unloading 
method, is shown to exhibit power law scaling with normal pressure, with the exponent of 
this relationship closely correlated to the fractal dimension of the surfaces. The electrical 
contact resistance at interfaces, measured using a controlled current method, revealed that 
the measured resistance is affected by testing current, mechanical loading, and surface 
topology. At a constant applied current, the electrical resistance as a function of applied 
normal stress is found to follow a power law within a certain range, the exponent of which 
is closely linked to surface topology. The correlation between stress-dependent electrical 
contact and normal contact stiffness is discussed based on simple scaling arguments. This 
study provides a first-order investigation connecting interfacial mechanical and electrical 
behaviour, applicable to studies of multiple components in energy systems. 

© 2016 Elsevier Ltd. All rights reserved. 


and adhesion [7-9]. In energy storage and conversion ap- 
plications the effective mechanical and electrical proper- 


Interfacial electro-mechanical behaviour is fundamen- 
tal to indicators of energy system performance such as 
electrical power loss, cycle life, and storage capacity in 
lithium-ion batteries [1,2], sodium-ion batteries [3], solid 
oxide fuel cells [4], photovoltaics [5] and thermoelectric 
systems [6]. Surface morphology plays an essential role in 
determining how contacting solids interact with one an- 
other in a variety of processes including thermal swelling, 
electrical conduction, electrochemical reactions, friction, 
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ties of granular electrode structures can be connected to 
microstructural characteristics [2,10]. However, the inter- 
facial properties in the existing modelling approaches are 
usually simplified [11]. 

Energy losses at interfaces are usually associated with 
ohmic heating (also known as Joule heating) due to the pas- 
sage of an electrical current through contacting surfaces. 
In the context of energy management, improved electrical 
contacts play a prominent role in mitigating energy losses 
in battery assemblies. The energy loss due to the electrical 
contact resistance (ECR) at interfaces between electrode 
layers and at contacts between electrodes and current- 
collectors can be as high as 20% of the total energy flow 
of the batteries under normal operating conditions [12,13]. 
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The effects of the mechanical properties and surface rough- 
ness of electrical contacts on the performance of electri- 
cal connectors are of great importance in terms of poten- 
tial drop and heat accumulation in contact zones [14,15].A 
significant increase in ECR can be caused by interfacial re- 
sistance due to the inevitable presence of resistive surface 
films, including corrosion deposits, fracture debris, oxide 
and hydrated layers at electrical contacts, resulting in ex- 
cessive ohmic heating. In extreme cases, the heat can bring 
about system failure through sparks, fire and even melting 
of system components |[12,16,17]. 

The stress dependence of ECR at rough surfaces can be 
associated with the varying true interfacial contact area 
during system operation. However, the direct quantitative 
evaluation of real interfacial contact area between bodies 
through either experimental measurements or numerical 
simulations remains highly challenging due to the com- 
plex multi-scale morphologies exhibited by rough surface 
structures [18-20]. Significant difficulties remain in re- 
lating interfacial electro-mechanical properties to surface 
structure descriptors. 

Employing electrical measurement, nanoscale mechan- 
ical testing and surface morphology characterisation, we 
investigated interfacial normal contact stiffness and elec- 
trical conduction behaviour at rough interfaces with ran- 
dom multiscale morphologies. First, we conducted contact 
stiffness measurements using flat-tipped diamond nanoin- 
dentation tests on a set of rough surfaces. Then, we ex- 
amined the evolution of electrical conduction with vary- 
ing compressive loads. Based on these results, discussions 
are extended to the relationship between electrical con- 
tact conductance and contact stiffness. This study demon- 
strates the importance of a multi-physics understanding of 
the origins of the electro-mechanical behaviour at inter- 
faces in order to improve the reliability and performance 
of electrical contacts in energy systems. 


2. Theoretical background 


Compared with the apparent or nominal contact area, 
the true area of contact at an interface is considerably 
smaller due to the existence of surface roughness. As 
shown in Fig. 1, when an electric current is conducted 
between two contacting solids, the restricted contact 
area, which depends on the size and spatial distribution 
of contacting asperities, causes additional constriction 
resistance (known as the electrical contact resistance, 
ECR) [21]. In addition to the constriction resistance 
resulting from the limited areas of true contact at an 
interface, ECR is also affected by the existence of resistive 
surface films, such as oxide layers [8]. Theories of ECR 
have since been further developed to include the effects 
of elastic-plastic deformation of the contacting asperities 
due to applied forces, multi-scale surface topography, size 
effects, and the contribution of insulating films between 
contacting bodies | 18,22-25]. 

Current flowing through rough interfaces is scattered 
across a large number of micro-contacts of various 
geometries, which are often assumed to be circular in 
theoretical treatments [19,26]. The constriction resistance 
due to the convergence and divergence of current flow 





Fig. 1. Schematics of electrical conduction through a rough interface 
exhibiting multi-scale surface features. 


at a single contact is represented in Fig. 1. The resistance 
of a single contact is dictated by the dominant electronic 
transport mechanism, which depends on the contact area 
and structure. When the radius of the micro-contact, r, is 
comparable or smaller than the average electron mean free 
path, A, the constriction resistance is dominated by the 
Sharvin mechanism, in which electrons travel ballistically 
across the micro-contacts. The resistance of a contact with 
area, a, is given by [27]. 


_ A(p1 + 2) (1) 
7 2a 


where o; and p2 are the specific resistivities of the 
contacting surfaces. On the other hand, whenr > 4A, 
the electron transport through the contact can be treated 
classically (Holm contact). The resistance can be expressed 
in the following form [8]: 


_ Jt (p1 + £2) (2) 
— Ae 


The total electrical conductivity, G,, of an interface is 
assumed to be the sum of individual conductivity G; = 
R, ' at micro-contacts, corresponding to the restriction 
resistances in parallel: 


G. = > G. (3) 


In the case of rough surfaces with multiple contacting 
asperities, there is a distribution in the size of the contact 
area. The Sharvin and Holm expressions should therefore 
be considered as limiting cases. 

In general, the area of contacts used in Eqs. (1) and 
(2), and therefore the contact resistance, depends on 
the applied pressure. Using theoretical and numerical 
approaches | 18,22,28,29], power-law type semi-empirical 
correlations between the contact resistance and the 
normal pressure have been proposed for rough interfaces. 
In particular, previous theoretical studies found the contact 
conductance to be linearly proportional to the incremental 
stiffness [18,22]. 

It should be noted that many mechanisms of surface 
Structure evolution have been observed during electri- 
cal conduction through rough interfaces, including dielec- 
tric breakdown of oxide layers, localised current-induced 
welding, chemical disorder arising with random compo- 
sition and oxidation processes in corrosive environments, 
and surface diffusion [30,31]. These phenomena are out- 
side the scope of this paper. 
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Fig. 2. SEM images and typical surface profiles of aluminium samples subjected to different surface treatments: (a) polished, S1; (b) sand blasted with 
300 «um-sized glass beads, S2; (c) sand blasted with 50 tum-sized glass beads, S3. 


Table 1 

Sample surface characteristics for different surface treatments. 
Sample RMS Fractal dimension, RMS slope 
type roughness/jum Dr 
S1 0.057 + 0.005 2.093 + 0.062 0.009 + 0.001 
S2 4.179 + 0.194 2.551 + 0.022 0.224 + 0.015 
S3 2.970 + 0.276 2.626 + 0.017 0.202 + 0.010 


3. Surface preparation and characterisation 


Round disks, with a diameter of 25 mm, made of 
aluminium alloy 5005 were used to fabricate specimens for 
both the measurement of the interfacial contact stiffness 
and ECR. For each individual sample, both the top and 
bottom surfaces were subjected to the same treatment 
using standard polishing and sand blasting procedures. 
The average diameters of the two selected groups of 
glass beads used in blasting treatments were 50 wm 
and 300 im. The sand blasting process was conducted 
for one minute, a duration which was sufficient to 
yield homogeneously and isotropically modified surface 
features. The sample surfaces were fabricated in such 
a way that each set of surfaces exhibited a distinct 
combination of surface roughness indicators, namely root 
mean squared (RMS) roughness and fractal dimension. 
Fig. 2 shows scanning electron microscope (SEM) images 
and typical surface profiles of the different surface types 
used in this work. Based on the three-dimensional digitised 
topographies obtained by optical surface profilometry 
(NanoMap 1000WLI), the mean values of RMS roughness, 
fractal dimension and RMS slope with standard deviations 
over ten scans on different samples were calculated, as 
shown in Table 1. These values are found to be comparable 
with descriptors of naturally occurring surfaces [32]. 
Values of RMS roughness were calculated as the RMS 
average of the profile height from the scanning. In 
the digitised scanning, the slopes of triangular units 
formed by three adjacent pixels are used to calculate 


Contact stiffness, 


Electrical conductance, Exponent ratio, a¢/Q 


QE AG 

0.463 + 0.022 0.816 + 0.081 1.762 + 0.194 
0.569 + 0.029 1.026 + 0.049 1.803 + 0.126 
0.605 + 0.022 1.494 + 0.134 2.469 + 0.239 


the RMS slope, which is commonly chosen as a higher 
order surface descriptor [33,34]. The scaled triangulation 
method [34] was used for the calculation of fractal 
dimension values. It was found that the smaller the 
particles used to modify the surfaces the larger the 
fractal dimension was. The fractal dimension, a cross- 
scale surface descriptor that incorporates localised and 
macroscopic surface information provides an effective 
means for modelling engineering surfaces with random 
self-affine multi-scale properties in the characterisation of 
surfaces and particles [35]. The advantage of using surface 
fractality as a cross-scale surface descriptor stems in part 
from the tendency of first order descriptors (e.g., maximum 
height or mean roughness of the surface) to be dominated 
by highest scale features, while secondary descriptors 
(e.g., slope) tend to be dominated by finest scale surface 
characteristics [33,34]. 


4. Contact stiffness at rough surfaces 


The surface contact stiffness of aluminium samples 
with different surface morphology was assessed using 
nanoindentation (Agilent G200) with three flat indenter 
tips of different diameters of 54.1 wm, 108.7 wm, 
and 502.6 wm (FLT-D050, FLT-D100, and FLT-D500, 
respectively, SYNTON-MDP, Switzerland). The reason for 
choosing flat tips is that the apparent contact area under 
the tip does not change with respect to the indentation 
depth, which is not the case for spherical or Berkovich 
tips. When the flat indenter tip first comes into contact 
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Fig. 3. Typical loading-displacement curves of nanoindentation tests on 
three types of surfaces. Ten partial unloading tests were carried out to 
isolate elastic contributions to contact stiffness under different loading 
levels. The three flat indenter tips used in the experiments (FLT-D050, 
FLT-D100, FLT-D500) are also illustrated for comparison. 


with the testing sample, the actual contact area is only a 
small fraction of the nominal contact area. The asperities 
of the sample surface at contact regions are then squeezed 
against the flat tip as indentation progresses as is shown 
in Fig. 3. In order to evaluate only the elastic responses, 
partial unloading tests were successively performed at 
ten intervals by decreasing the applied load by 10% each 
time. The loading level of each subsequent unloading stage 
is twice that of the previous unloading stage, with a 
maximum load of 500 mN during the last unloading step. 

Mean stiffness values were obtained by averaging data 
of ten indentation tests at different locations for each 
surface type. The unloading stiffness is here defined as the 
initial slope of the unloading curve, k = dF/dS, where F 
designates the normal force and S is the indentation depth. 
Subsequently, the reduced elastic modulus E, was derived 
from the measured unloading stiffness as 


2 

k BEA, (4) 
where A is the apparent contact area of the indenter 
tip and #6 is a geometrical constant, taken as unity for 
a flat punch [36]. Eq. (4) is a fundamental equation for 
assessing the elastic properties in nanoindentation tests. 
The reduced modulus depends on the elastic properties of 
both the tested specimen and the indenter tip: 


1 1-vue 1-vu 


E, E. E; 
where vu; and vu, represent the Poisson’s ratios of the 
indenter tip material and the tested specimen respectively. 
For the diamond indenter tips used in this research, E; and 
vu; are typically 1140 GPa and 0.07, respectively. Eqs. (4) 
and (5) allow the estimation of E, from measured values 
of A and k, while for vu. we simply use the Poisson’s ratio 
v* of bulk aluminium (uv, = v* = 0.3). 

By using different sized flat tips, the stress range 
extends over several orders of magnitude. With the same 
maximum force (500 mN) provided by the nanoindenter, 
the maximum stress produced with FLT-D050 was around 
100 times larger than that produced with FLT-D500. The 


; (5) 








Stress provided by all the three indenter tips ranged from 
0.005 MPa to 200 MPa, spanning five orders of magnitude. 
The contact stiffness measured over this range of applied 
stresses varies approximately from 0.01 GPa to 55 GPa. 

Fig. 4 shows the evolution of contact stiffness with ap- 
plied force for the different surfaces. Here, we normalised 
the contact elastic modulus E, by the Young’s modulus of 
aluminium alloy 5005, E* = 69.5 GPa. The force is nor- 
malised by E*A, where A is the projected area of the cor- 
responding tip. The measured contact stiffness increases 
with the loading force, for all tested samples. At the same 
applied stress level, the surfaces after sand blasting treat- 
ment (samples 2 and 3) show a smaller value of contact 
stiffness with respect to that of the polished surface (sam- 
ple 1). The surface blasted with glass beads of 50 wm di- 
ameter (sample 3) presents the lowest contact stiffness of 
all the three types of surfaces. 

We express here the power-law relation of the contact 
stiffness with the applied normal force 


E, x (F)"*, (6) 


where ag is the exponent of the power-law function 
[37,38]. It should be noted that the fitting curves are 
achieved excluding the contributions from the measured 
stiffness under stress levels higher than 100 MPa, where 
the surface shows an apparent yield phenomenon. For all 
the three surface types, the value of the exponent a; varies 
from 0.4626 to 0.6048 (in Table 1), changing as the fractal 
dimension increases. In comparison, the typical value in 
cases of Hertzian contact of two elastic spheres is 1/3, 
as shown in Section 6. The power-law relationship found 
here experimentally is in good agreement with previous 
theoretical predictions on a quantitative basis [18,37,38]. 


5. Electrical conductance at rough surfaces 


For each surface type, interfacial electrical conductance 
was measured for stacks of eleven disks, giving ten rough- 
to-rough interfaces. Analysis was achieved by means of a 
source/measurement unit (SMU B2900A, Agilent) across 
a range of applied compressive loading forces. In this 
experimental setting, we measured the resistance created 
by ten interfaces instead of a single interface, aiming to 
achieve a higher precision, larger linear range and better 
robustness against the measurement noises from the 
connecting wires, loading device and measurement unit. 
Using multiple interfaces further reduces experimental 
errors arising from inhomogeneity in surface treatment 
processes. 

Prior to the measurement of force-dependent resis- 
tance, we performed resistance creep tests and sweeping 
current tests to select the most appropriate testing current 
and time to minimise influences on the measurement from 
the applied current. Full procedures and results have been 
previously published in greater detail [17]. The applied 
Sweep current test consisted of two phases: a “loading” 
phase (P1) with current increasing logarithmically from 
0.0001 A to 1.5 A, followed by an “unloading” phase (P2) 
with current decreasing logarithmically from 1.5 A back to 
0.0001 A. Both phases were conducted under conditions of 
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Fig. 4. Curve fitting for the normalised stiffness, E,/E*, and the normalised applied force, F/(E*A), for three tested surfaces, with E* being the Young’s 


modulus of the tested material, and A the apparent contact area. 
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Fig. 5. Typical measured results for polished samples using current sweep under various stresses (0.031 MPa, 0.061 MPa, 0.122 MPa, 0.245 MPa and 0.490 
MPa, corresponding to loops 1-5, respectively) with solid lines representing the first phase (P1) and the dashed lines showing the second phase (P2). 


constant normal load. During the sweeping loops, the volt- 
age was recorded at a frequency of 2 kHz. The two-phase 
Sweeping process was completed within 0.2 s in order to 
avoid significant time dependent resistance degradation. 
Fig. 5 shows the typical resistance-current characteris- 
tics for polished samples obtained from sweeping current 
tests. Each individual loop corresponds to a distinct load. 
The five loops shown demonstrate similar trends known as 
the Branly effect [30,31], i.e., the measured resistance be- 
gins to drop irreversibly after the testing current reaches 
a certain value. The process is featured by voltage creep, 
hysteresis loops, and voltage saturation effects [31,39]. The 
corresponding threshold current values for loops (1-3) are 
approximately 150 mA, 200 mA and 400 mA, respectively, 
and the value seems to be positively correlated with the ap- 
plied normal load. However, the Branly effect tends to be 
harder to capture at sufficiently high stress levels, shown 
in loops (4-5). For all five loops, when the testing cur- 
rent is higher than approximately 5 mA and lower than 


the threshold current values, the measured resistances re- 
main stable at two plateaus in both P1 and P2, and can 
therefore be defined as ohmic resistance (the testing cur- 
rent is directly proportional to the measured voltage). At 
low testing currents (lower than 1 mA), the measured re- 
sistance exhibits instability with the prominent measure- 
ment noises. The measured resistance obtained from sub- 
sequent sweeping tests will follow the path of the unload- 
ing phase (shown in the dashed lines in Fig. 5) [17]. 

The experimental results in Fig. 5 indicate that both 
mechanical loading and electrical current alter the surface 
morphology and broaden the gap in measured resistance 
between P1 and P2. The contacting asperities can be 
regarded as a resistor network changing with the applied 
current, mechanical load, and measurement time. At the 
interfaces, the electrical current results in the physical and 
chemical modification of sample surfaces, which involves 
many processes, including the rupture of the oxide layer 
due to compression, and the localised heating induced 
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Fig. 6. Stress-dependent electrical conductance of different surfaces under various loading levels with a testing current of 10 mA, with E* being the value 
of the Young’s Modulus of the tested material, A the projected area of the tested samples, and Gp being set to 1 2~'. The curve fitting were conducted for 


loading levels in a range of [0.031 MPa, 3.973 MPa]. 


by current. A high level of applied stress leads to better 
stability and repeatability of ECR measurements [17,31]. 

Based on the performed sweeping current tests the 
electrical resistance was measured under various stresses, 
for samples exhibiting different surface morphologies. For 
each type of sample, five series of tests were conducted and 
the resistances were evaluated at 16 different stress lev- 
els from 0.020 MPa to 8.936 MPa. The measured time was 
0.01 s for each individual data point, in order to avoid sig- 
nificant effects arising from ohmic heating and associated 
time dependent resistance degradation. The testing cur- 
rent was set at 10 mA, where all the three types of samples 
display an ohmic behaviour under varying electrical and 
mechanical loads. The interfacial electrical contact conduc- 
tance was subsequently calculated through G, = 1/(R, — 
Ro), where Ro (~0.06 &2) is the combined resistance of the 
bulk material of identical size as the disk stack (~2.53 Q), 
wires and connections used in the experimental setting. 

As shown in Fig. 6, the measured conductance of disk 
Stacks increases considerably with pressure, converging 
to a value close to the bulk conductance of the material. 
For given stress levels (<4 MPa), samples blasted with 
50 wm sized glass beads (S3) usually present the lowest 
conductance among all the three types of samples. At low 
levels of applied stress (less than 0.5 MPa), the conductance 
is spread across a wider range. Similar to Eq. (6), we use 
a power-law function to express the correlation of the 
conductance with applied normal load as 


Ge a (F)"S (7) 


By fitting the conductance/pressure curves from 0.031 MPa 
to 3.973 MPa, the power law exponent az is found to be 
0.816, 1.026 and 1.494 respectively for polished surfaces 
(S1), surface blasted with 300 wm particles (S2) and those 
treated with 50 wm particles (S3). The exponent values 
increase with the fractal dimension, shown in Table 1. 
Moreover, for all the three types of surfaces as shown 
in Fig. 6, the electrical conductance reaches a plateau 


under higher stresses, with the plateau value correlating 
to the RMS roughness. In the lower stress regime, the 
experimental data no longer seems to follow the power 
law. 


6. Discussion 


The key experimental results for contact stiffness 
and electrical conductance, measured for three types of 
rough surfaces, are summarised in Table 1. Both the 
contact stiffness and electrical conductance increase with 
the applied force, exhibiting power law behaviours with 
exponents a, and dc, respectively. These exponents vary 
with the surface roughness and increase with the fractal 
dimension. In contrast, no evident correlation between the 
RMS roughness value and the exponent was found. This 
Suggests that the correlation between contact stiffness, 
electrical conductance and applied force is dominated by 
fine scale surface characteristics. 

We rationalise the experimental findings by developing 
the following scaling arguments. Both the contact stiffness 
and conductance primarily depend on the true contact area 
A,-, which evolves during mechanical loading and cannot 
be determined in a direct way based on the considered 
measurement methods. As a workaround, we estimate the 
true contact area based on the following expression for the 
incremental stiffness: 


2 
k= pak Ao (8) 


where E” is the (constant) reduced elastic modulus calcu- 


lated for the bulk elastic properties of the tested mate- 


rial and indenter: E’, = ((1 — v**) /E* + (1 — v7) /Ei) 


and f’ is a geometric factor of the order of unity. By writ- 
ing Eq. (8), we assume that the effect of surface roughness 
on the measured incremental stiffness can be described 
by considering the true contact area A, (rather than the 
project contact area A) and bulk material properties in the 
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fundamental Eq. (4). By comparing Eq. (4) and (8), and con- 
sidering that E; >> E* > E,, we obtain the following scaling 
relation: 


E./E* x B’\/A-/A. (9) 
Next, we consider the true contact area to be the sum 
of n individual contact areas, with average a = A,/n. 


Here, individual asperities are assumed not to interact 
with one another during deformation. In order to relate 
the evolution of a to the applied force, we first rely on a 
classical result of Hertzian contact theory. Representing a 
single contact by two spheres with radii R; and Rz squeezed 
against each other, the contact area varies with the applied 
force according to 


3RF \2/3 
C= 7 & ; (10) 
r 





where R = (1/R; + 1/R2)~' is the equivalent radius of the 
two spheres, and the reduced modulus E’ was introduced 
in Eq. (8). Eqs. (9) and (10) indicate that the contact 
stiffness is a power function of the load, with an exponent 
1/3. This simple scaling analysis is not consistent with 
our experimental findings for w;, which takes significantly 
higher values. However, the scaling analysis based on the 
Hertzian contact theory does not consider the changing 
number of contact asperities, n, for the increasing load. 
Furthermore, at the rough interface, the contact areas are 
not uniformly distributed [40], and interactions between 
asperities can exhibit complex deformation mechanisms, 
such as plastic deformation, adhesion, and friction. 

On the other hand, introducing relation (10) into 
Eqs. (1) and (2) for the Holm and Sharvin resistance at a 
single contact, one finds that: 


4 (3RF\*° Qn (3RF\2° 
— \4E’, Ap \ 4E’, 


where 90 = (0; + (2). Combining (11) with (3), we 
find that the total conductance of the rough surface, G,, 
approximately scales with the force following a power 
law with the exponent ranging from 1/3 (Holm) to 
2/3 (Sharvin), depending on the dominant conduction 
mechanisms at individual contacts. 

We further consider the contact model for a conical 
punch [41], where the contact area, a, is found to be linearly 
proportional to the applied force, F. With the same analysis 
as above, an exponent a; = 0.5 can be derived for the con- 
tact stiffness, and an exponent a, = 0.5 ~ 1 for the elec- 
trical conductance. This provides a better representation 
for the exponents of contact stiffness and electrical con- 
ductance as compared to the prediction by the Hertzian so- 
lution. Again, the exponents derived from this simple scal- 
ing analysis are lower than the experimental values of ac. 
This may also be due to the fact that the scaling neglects 
the increase in number of contacting points under increas- 
ing compression. 

Despite these discrepancies, it is interesting to consider 
the ratio of the exponents for contact stiffness and 
electrical conductance, a@¢/a,. This ratio characterises the 
power law relation between the conductance and contact 








stiffness, with G. « (E-)%¢/“£. According to the scaling 
analysis, this ratio ranges from 1 (Holm mechanism) to 
2 (Sharvin mechanism). Experimentally, an approximate 
value of 2 was found for loading levels in a range of 
F/(E*A) € [5 x 10-’,5 x 10°]. For sample 3, similar 
fitting in the low load region gives a higher value of ac, 
and hence a higher value of the ratio a¢/ag, suggesting 
a larger proportion of Sharvin-type contacts at low loads. 
Similar observations seem to hold as well for samples 1 
and 2, but the transition takes place at even lower loads. 
As the load increases, new asperities come into contact, the 
contacting points enlarge and small microcontacts merge 
forming large contacts, resulting in better conduction. 
The dominant conduction mechanism transitions from 
a Sharvin-type to a Holm-type with the exponent ratio 
decreasing from 2 to 1. Under sufficiently high forces, 
and hence high contact areas, the electric and mechanical 
properties converge to those of the bulk material, as 
expected. 

The ratio w/a, also tends to increase with the fractal 
dimension. A surface with a higher fractal dimension 
demonstrates Sharvin dominated conductance (a@¢/ag ~ 
2), while a less fractal surface presents combined Sharvin 
and Holm-type conductance (a@¢/ag between 1 and 2). 

Note that in the contact stiffness measurements, 
the flat indenter tips can be considered as rigid flat 
surfaces (E; >> E*), corresponding to a rough-to-flat 
contact problem. In comparison, our interfacial electrical 
resistance experiments involve rough-to-rough contacts. 
However, a scaling analysis based on rough-to-flat contact 
would yield identical exponents in the power law functions 
(10) and (11) [28,33]. 

In the experiments, both contact stiffness and conduc- 
tance may be affected by oxide layers at the sample sur- 
faces. Aluminium alloys ubiquitously exhibit thin passi- 
vated hydrous and oxide layers arising from reaction with 
atmospheric oxygen and water. This nanoscale layer ex- 
hibits locally divergent mechanical properties in a region of 
thickness typically less than 10 nm, which is significantly 
less than the depth of indentation performed in the cur- 
rent work. The influence of oxide layers is thus expected to 
be of limited significance in the present contact mechan- 
ics study. In the analysis of ECR behaviour, the oxide layer 
acts as an insulator. However, due to its limited thickness, 
the measured conductance is only sensitive to the presence 
of this layer at lower loads. For this reason large measure- 
ment uncertainties are evident at low loads with the mag- 
nitude of these fluctuations dependent on specimen sur- 
face structure, as shown in Fig. 6. Therefore in this study 
the effect of the oxide layer is minimal and does not inter- 
fere with the findings. 

The observations made in this study can provide in- 
sights into the physical origin of the topological depen- 
dence of transport phenomena in energy materials applied 
in conversion, storage and generation systems. Parametric 
Studies in to the performance of energy systems often yield 
unexpected behaviour arising from changes to the struc- 
ture or processing of complex materials such as granular 
electrodes [2,10]. The present work suggests that the struc- 
ture and mechanics of interfaces in these systems may be 
in part a contributing factor to the observed processing de- 
pendence of performance. 
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7. Conclusion 


We performed experimental investigations into the 
contact stiffness and electrical contact resistance at rough 
interfaces, with a specific focus on their dependence on 
applied force. The change of these interfacial electro- 
mechanical properties under different loading conditions 
can be associated with changes in the true area of 
interfacial contact. The measured contact stiffness and 
electrical conductance have been found to exhibit power 
law relationships with normal pressure across a wide 
range of applied stress, expressed as E- «x (Fy)“- and 
G. « (Fy)*¢, respectively. The corresponding exponents 
of these relationships were found to be closely correlated 
to surface fractality with the absolute values of a; and 
Qc, increasing with the fractal dimension of the surfaces. 
The presented experiments on load-dependent contact 
stiffness and electrical contact resistance provide an initial 
step towards connecting interfacial electro-mechanical 
properties and surface topology, which is of value in 
interpreting the properties of various energy materials 
and components. Further investigation is warranted to 
fully understand these phenomena and interpret the 
interface-morphological dependence of energy material 
performance. 
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5 Electrical Transport in Binary 
Percolation Networks 


Granular materials are actually multi-contact systems with particles. The electrical transport in many 
types of granular materials is realized by current paths formed by the inter-particle contact, the 
electrical properties of which depends on the inter-particle force and the interface structure, as shown 
in Chapter 3 and 4. The overall conduction behaviour can be simulated by resistor-capacitor networks. 
Specifically, each individual network node is a particle of the considered granular packing, and each 
individual network element, such as a capacitor or a resistor, represents a contact in the granular 
materials. In this chapter containing Paper 5, we studied the influences of network geometry and 
boundary conditions on the electrical responses of random binary percolation networks. For varying 
network configurations in terms of network aspect ratio and electrode dimension, a unified description 
focusing on the centre and the span of the emergent scaling region is provided for the universal 
dielectric responses. The findings of this paper can guide the design and testing of finite disordered 
systems with multi-contacts. 


This paper has been published in PLoS ONE in 2017. I was the primary researcher and author of this 
paper, being supervised by Dr. Yixiang Gan, Dr. Dorian Hanaor. 
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Abstract 


In this paper we apply lattice models of finite binary percolation networks to examine the 
effects of network configuration on macroscopic network responses. We consider both 
square and rectangular lattice structures in which bonds between nodes are randomly 
assigned to be either resistors or capacitors. Results show that for given network geome- 
tries, the overall normalised frequency-dependent electrical conductivities for different 
capacitor proportions are found to converge at a characteristic frequency. Networks with 
sufficiently large size tend to share the same convergence point uninfluenced by the bound- 
ary and electrode conditions, can be then regarded as homogeneous media. For these net- 
works, the span of the emergent scaling region is found to be primarily determined by the 
smaller network dimension (width or length). This study identifies the applicability of power- 
law scaling in random two phase systems of different topological configurations. This under- 
standing has implications in the design and testing of disordered systems in diverse 
applications. 


Introduction 


The bulk behavior of complex systems comprising disordered multi-phase components is of 
importance in diverse applications including supercapacitors and batteries [1-4], dielectric 
material characterization [5-10], the mechanics of structures [11-13], fracturing process [14], 
thermal analysis [15] and soil probing [16]. In such systems, various parameters govern the 
electrical, thermal, chemical, and/or mechanical properties of components of a system across 
multiple scales from molecular up to macroscopic length-scale. Experimental and computa- 
tional research efforts are increasingly conducted in order to gain insights into the manner in 
which these properties combine across scales to determine overall system performance. 

In particular, the AC conductivity of systems that can be schematically represented as mix- 
tures of electrical components has been the subject of numerous investigations that have 
shown power-law scaling with frequency arising through different relaxation mechanisms [11, 
17-19]. Above a critical frequency this scaling of AC conductivity is described by Jonscher’s 
power law[20] and has been experimentally observed across diverse conductor-dielectric 
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composites and porous materials [5, 7, 18, 21, 22] with this scaling being termed the “Universal 
Dielectric Response” [17, 18, 21, 23]. This emergent property does not arise directly from any 
particular physical or chemical properties of the involved components, but rather is a conse- 
quence of the way components combine [11, 19, 21, 24, 25]. Such dielectric mixtures have 
been effectively approximated as a random network of resistors and capacitors [18, 19, 24] 
with representative conductors exhibiting a constant conductance 1/R and dielectric compo- 
nents exhibiting a variable complex admittance iwC, which is directily proportional to an 
angular frequency w, as illustrated in Fig 1. Useful asymptotic formula for the emergent net- 
work admittance including both the effects of component proportions and the network size 
can be obtained based on the spectral method [21, 26] and the averaging approach [11, 21, 27]. 
However, establishing a more rigorous estimation necessitates numerical analysis. 

From previous numerical studies [11, 17-19, 21, 23, 25, 28], the typical obtained conductiv- 
ity-frequency spectrum of a square lattice resistor-capacitor (RC) network can be divided into 
three regions of angular frequency, w, governed by the proportion of capacitors, p,, and net- 
work size, N (a) an emergent region for intermediate frequencies; (b) two percolation regions; 
(c) transitions between the two above-mentioned [21]. Symmetry is found between the overall 
responses at high and low frequencies, which can be correlated to percolation behavior [26, 29, 
30]. For low and high frequencies, current tends to percolate predominantly through resistors 
and capacitors respectively, as these will exhibit relatively lower impedance, as presented in Fig 
1. For intermediate frequencies, where the values of admittance for the resistors and the 





L 


Fig 1. A lattice network containing Wx L =5 x 5 randomly distributed resistors and capacitors. The network width Wand 
length L indicate the numbers of horizontal and vertical elements in a single chain, respectively. The values of By = 2, B; = 4, and 
Bo = 6 present the number of nodes connected directly to the electrodes. Percolation paths formed by resistors and capacitors are 
shown in thick blue and orange lines, corresponding respectively to the dominating modes at low and high frequencies, for the Bo 
configuration. 


doi:10.1371/journal.pone.0172298.g001 
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capacitors are close, we observe the power-law emergent behavior whereby conductivity is 
proportional to w®%, with a x p [11, 21, 27]. 

Motivated by the applicability of these networks for representing real-world disordered sys- 
tems, in this paper we reexamine the universality of the emergent power-law scaling observed 
in previous work by further considering the significance of the network aspect ratio and 
boundary conditions on the convergence point and the span of the emergent region, the two 
key network characteristics of universal scaling behavior. 


Methods 


In this paper, we extend the square lattice RC networks to rectangular ones with N= W x L ele- 
ments distributed between two bus-bars, one of which is grounded and the other raised to a 
potential, | Vole’. A system of complex number linear equations is set up by applying Kirchh- 
off’s current law (for complex currents) on each individual node of the RC network. For the 


node k, we get 


I,(t) _ ZT, [eortew = =? (Vielen _ Vy let Zo = 0, (1) 


where J; is the sum of the currents (negative or possitive) flowing towards the node k, from 
connected components. The impedance of a component connected to node k, Z; ;, is random- 
ized to be either R or 1/iwC. The voltage potential of the connected node j, Viger, with 
respect to that of the node k, is represented as complex-valued function of time, ¢, with 9, ; 
being the relative phase. The value of n is determined by the location in the network, equaling 
the number of connected components. More specifically, n = 4 for ones located away from the 
boundaries, n = 2 for the lattice corners, n = 3 for the nodes on the boundaries excluding cor- 
ners. The two electrodes are also regarded as nodes with n = B. The electrode dimension, B, is 
defined as the number of nodes connected directly to the electrodes, e.g., B= W + 1 when all 
the elements along the boundary side are connected to the electrode. Each single node is repre- 
sented by a corresponding linear algebraic equation, resulting in (W + 1) x (L + 1) equations 
for all the nodes. Two additional equations can be obtained from the electrodes. By solving 
these equations, the potential of each node and the current going through each bonds in the 
network can be calculated. Thus, for a given applied potential difference between electrodes, 
| Vole“, the macroscopic admittance can be given by Y = |J,|e't# /|V, |e", where |, |e*+#” 
is the obtained overall current flowing into the ground. Additionally, Frank-Lobb techniques 
can be employed to reduce the network size, thus improving the computational efficiency [31]. 

Here, the network aspect ratio and the electrode dimension are considered as variables, in 
order to investigate the influence of network configuration and boundary conditions on mac- 
roscopic responses. The network size is determined by its width, W, and length, L, which rep- 
resent the numbers of components in a single chain along the horizontal and vertical 
directions, respectively. In the network circuit, two electrodes of identical dimension are con- 
nected to elements located symmetrically in the center of the two vertical boundaries. The fre- 
quency-dependent macroscopic responses obtained from different configurations, in terms of 
network length, L, and width, W, are normalized through 

\Y|JRL 


Y= , © =aRC, (2) 
W-+1 





where R and C are the resistance and capacitance values of resistors and capacitors in the con- 
sidered network, respectively. This normalization process is applied in order to include the sig- 
nificance of all the elements in the overall network behavior, represented by the equivalent 
admittance, by considering rules for simple series and parallel combinations of components. 
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Fig 2. Normalised admittance module as a function of frequency. Numerical results obtained from three groups of differently 
configured networks (denoted as W~ L_B) are presented, with capacitor proportions, p,, varying from 0% (corresponding curves are 
shown in black) to, 25% (red), 50% (blue), 75% (green), and to 100% (brown). The phase responses are depicted in the inset. For 
each network configuration with a given capacitor proportion, five simulations have been realized. 


doi:10.1371/journal.pone.0172298.g002 


The span of the emergent region, S, is defined as the horizontal distance between the intersec- 
tions of the power-law function, y = w”, with top and bottom percolation admittances, Y, 


and Y, (averaged from multiple simulations), for low and high frequencies, respectively, as 
shown in Fig 2. 

The network behaviors shown in Fig 2 in the frequency domain are primarily governed by 
percolation effects [32], which are closely linked to the frequency-dependent conductivity of 
each single bond in the network. In the studied rectangular RC networks with two types of 
bonds (i.e., resistor and capacitor, the admittance ratio of capacitor elements with respect to 
resistors is iwRC) have been considered to describe the responses of random binary networks. 
The observed universal scaling behavior in Fig 2 can be also found in other networks 
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containing two types of elements, indexed by a and J, exhibiting differences which can be 
described in the form of S, = w*S;, e.g., mechanical stiffness, thermal conductivity, and chemi- 
cal reaction rate, etc. [11, 14, 25]. 


Results and analysis 


The AC electrical responses of three groups of different sized networks (expressed in the form 
of W x L_B, representing network widthxlength_electrode dimension: 20 x 20_21, 100 x 
20_101, and 20 x 100_21) with varying p, are plotted in Fig 2. The three types of regions 
observed in the conductance spectroscopy of square networks can be also observed here in 
rectangular ones. For a given network geometry, the obtained admittance spectroscopies for 
various p,, are found to intersect at a convergence point with the characteristic frequency of 

w = 1/RC where resistors and capacitors contribute equally to the overall conduction. This 
point also appears to be the center of the emergent region. The normalized characteristic 
admittance (the values of Y at the characteristic frequency) at the convergence point is close to 
1. This indicates that the network at the characteristic frequency perform effectively as a 
mono-element network, which has a phase angle reaching the extremum value, as is shown in 
Fig 2. It is further evidenced that the normalized emergent regions of the three groups of net- 
works coincide with each other presenting universal features. The normalized admittance in 
this common emergent region appears to be uninfluenced by the network aspect ratio (length/ 
width) or the electrode dimension. However, differences can be found at percolation regions 
along with the corresponding transition regions. Variation of the width or length can poten- 
tially change the percolation thresholds which will determine the responses, following a resis- 
tive-percolated (plateaued) or capacitive-percolated (upwards or downwards) trend, 
corresponding to the low and high frequency ranges, respectively. 

Statistical analysis of the normalised characteristic admittance for different-sized square 
networks (from 5 x 5 to 600 x 600) with B = W + 1 was conducted and the standard deviation 
(STD) of the normalized characteristic admittance is presented in Fig 3. It is found that all val- 
ues of normalized characteristic admittance obtained with various p, (averaging over ten simu- 
lations) are in the range of (0.95, 1.05) for networks with more than 10 x 10 elements. The 
variance tends to diminish as the network size increases. This can be explained by considering 
boundary effects that relatively smaller networks have higher percentage of boundary elements 
(connected to four other elements rather than six in the bulk region). Responses of larger net- 
works perform with little influence from the boundary. For a given sized square network a 
larger variation is found for cases of p = 1/2, as such conditions lead to an equal likelihood of 
resistive-percolated and capacitive-percolated network responses at low and high frequencies. 
Consequently, there are four possible qualitatively different types of response for any realiza- 
tion of the system.[21] Different available responses potentially introduce dispersion and 
uncertainty of the network behavior in both percolation and emergent regions, as can be seen 
in Fig 3. 

We consider 2D rectangular networks with various L/W and B/(W + 1) ratios, in order to 
study the effects of network size and electrode dimension on frequency dependent responses. 
The variations of normalized characteristic admittance for rectangular networks (not shown) 
are comparable to those of square networks. Here, the convergence-divergence behavior is 
tested with three groups of rectangular networks, which have the fixed width, W, of 20, 50, or 
100 elements, respectively. The results obtained from the three groups coincide with each 
other, and typical results for L/W = 1.5, 0.8, 0.2 are shown in the Fig 4A. The contour of nor- 
malized characteristic admittance shown in Fig 4B presents a clear trend approaching 1, as the 
network length and the electrode dimension increase. For an RC network with a given size, a 
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Fig 3. Mapping of the standard deviation of normalised characteristic admittance values. For varying capacitor proportions 
from 0.1 to 1.0, different-sized square networks (from 5 x 5 to 600 x 600) were considered. For a given network size and capacitor 
proportion, ten RC networks were generated and used in the simulations to obtain the averaged normalised characteristic 
admittance, represented by the black dot. The STD of these points are used for mapping with the colour indicating the STD values, as 
detailed in the legend. 
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smaller value of electrode dimension tends to constrict the current to fewer paths at zones near 
the electrodes, thus, effectively reducing the network length. However, this influence will 
diminish as the network lengthens. Networks with large enough length tend to share the same 
intersection point uninfluenced by the boundary condition. In this case, these networks can be 
defined as homogeneous systems, with the responses in emergent region unaffected by the net- 
work configuration and the electrode dimension. A smaller length/width ratio or electrode 
dimension will usually lead to normalized characteristic admittance values smaller than 1. 
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Fig 4. Influences of network aspect ratio and electrode dimension on the values of characteristic 
admittance. (A) The dependence of characteristic admittance on the B/(W + 1) ratio, for cases of different 
widths, W, of 20, 50, and 100 but same L/W, including 0.2, 0.8, and 1.5 (corresponding to A, Mf, and @, 
respectively). (B) Mapping of characteristic admittance values obtained from different-sized networks with 
varying L/Wand B/(W + 1) ratios. The square network is marked by the black dot. 


doi:10.1371/journal.pone.0172298.g004 
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For cases of electrodes of various geometries on the boundary or embedded in RC net- 
works, the influence of the electrodes is mainly induced by the elements connected directly to 
the electrodes. As the network size increases, boundary elements and electrode-affected ele- 
ments will have a decreasing percentage. Therefore, the electrode size effect along with the 
boundary effect will be unobservable for a sufficiently large network. The trend can be found 
in Fig 4B that the characteristic admittance values for an increasingly large network locate well 
upwards and to the right from the red zone, asymptotically approaching 1. This universal char- 
acteristic can be also extended for an infinite network, with the normalized intersection value 
to be 1. 

The power law dependence of the electrical responses on frequency is of a universal nature 
for a wide range of complex RC networks. To further interpret this emergent behavior, we 
investigate the spans of the power-law emergent region. The convergence point tends to be the 
geometric center of all the emergent regions for various p,. By considering the region center 
reported here combined with the span, S, the emergent scaling behavior region can be well 
described. 

It has been observed [21] that, for a square network with p, = 0.5, the span of the power-law 
emergent region increases without bound as the network size increases. In this paper, using 
the results from square networks as references, we compare the spans of emergent regions 
obtained for various sized rectangular networks with different electrode dimensions, as p, var- 
ies. We found that when the normalized characteristic admittance value is sufficiently close to 
1, it is the smaller dimension, S,,;,, = min(W, L) that determines the span of emergent region, 
while the electrode dimension has little influence on the length. Different sized homogeneous 
networks with the same s,,,;, tend to effectively present identical responses for the emergent 
region with a given p,. However, discrepancies can be found for responses at low and high fre- 
quencies dominated by percolation behavior. This likely to be the case also for an intersection 
value far away from 1, but with lower accuracy due to the instability of the responses on 
account of the boundary effects. Only the results with p, from 0 to 0.5 are discussed and pre- 
sented in Fig 5, with the consideration of symmetricity of network responses. 

The results shown here for varied network and electrode dimensions shed light on the 
behavior of infinitely-large binary percolation-type network and large networks with irregular 
boundaries (e.g., in the shape of spline curves) and electrodes (e.g., with the geometry of circu- 
lar zones, embedded in the network, or unequal-sized electrodes). As long as the network size 
is significantly larger than the electrode and boundary dimension, the presenting universal fea- 
tures will not be influenced by the boundary and electrode conditions. This enables the net- 
work responses to reach a robust and reliable status at the emergent region with the span being 
determined by the effective network size in the order of D,” (D, is the distance between positive 
electrode and ground, indicating the shortest current path) and p, (dominating the slope of 
universal power law). Additionally, evident trend presented in Fig 5 supports that infinitely 
large emergent scaling regions can be observed for various capacitor proportions. 


Conclusion 


We studied the influences of network geometry and electrode dimension on the electrical 
responses of rectangular RC networks. The universal scaling behavior can be fully character- 
ized using the center and the length of the emergent region, i.e., the convergence point, and 
the span, respectively. For both square and rectangular networks, a convergence point is 
observed at the characteristic frequency, w = 1/RC, which usually appears to be the center of 
the emergent region. At this characteristic frequency, the normalized admittance value | Y|RL/ 
(W + 1) approaches 1 as the length-to-width ratio and electrode dimension increase. For a 
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Fig 5. Dependence of the emergent region span on network size and electrode dimension. The results for square networks are 
shown by solid lines with error bars obtained across ten simulations for each point. Results of rectangular networks (10 x 100_11, 
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Pc = 0.5 with those of rectangular networks with various network configurations. 


doi:10.1371/journal.pone.0172298.g005 





defined homogeneous network, the span of the emergent range is primarily determined by the 
shorter dimension of width and length. These observations provide a unified description for 
the emergent scaling properties of network responses for random two-phase systems with 
varying topological configurations. The comprehensive understanding of this emergent scaling 
can guide the design and testing of disordered systems in terms of determining testing condi- 
tions (e.g., the shape, size, location, and spacing of the fixtures), boundary conditions, and sys- 
tem dimensions. 
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6 Electrical Responses of Granular 
Materials 


In this chapter, we extend the studies of electrical conductance in granular materials under DC 
scenario to AC scenario over a wide range of frequencies. This chapter containing Manuscript 6 
presents experimentally observed universal AC scaling behaviour in conductive granular systems 
under different stress conditions. The spectra of impedance moduli on frequency follow a unique 
master curve with characteristic frequencies corresponding to the onset of conductance dispersion and 
measured direct-current resistance as scaling parameters demonstrating — stress-independent 
universality. Three-dimensional RC network topographies are extracted from the random granular 
packing realized by DEM simulation. The information collected at inter-particle contacts are scaled 
up to simulate the macroscopic responses of the granular packing through the formed network 
structures. The inherent electrical responses are found to be governed by the network topology and 
inter-particle contact morphologies formed under a given stress level, by means of introduced 
numerical model based on 3D packing structures and equivalent RC networks. In Chapter 6, we break 
the electrical transport in granular materials down into electrical transport at individual inter-particle 
contact, in order to utilise the knowledge of contact mechanics presented in Chapter 3 and 4. 
Subsequently, based on the studies shown in Chapter 5, the information at inter-particle contact level 
is integrated upwards to macroscale through 2D/3D regular/irregular RC network structures to reveal 
the physical origin of the AC responses of granular materials. 


I was the primary researcher and author of this manuscript, being supervised by Dr. Yixiang Gan, Dr. 
Dorian Hanaor. 
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Abstract 


We experimentally and numerically examine stress-dependent electrical transport in granular materials to 
elucidate the origins of their universal dielectric response. The ac responses of granular systems under 
varied compressive loadings consistently exhibit a transition from a resistive plateau at low frequencies to 
a state of nearly constant loss at high frequencies. By using characteristic frequencies corresponding to 
the onset of conductance dispersion and measured direct-current resistance as scaling parameters to 
normalize the measured impedance, results of the spectra under different stress states collapse onto a 
single master curve, revealing well-defined stress-independent universality. In order to model this 
electrical transport, a contact network 1s constructed on the basis of prescribed packing structures, which 
is then used to establish a resistor-capacitor network by considering interactions between individual 
particles. In this model the frequency-dependent network response meaningfully reproduces the 
experimentally observed master curve exhibited by granular materials under various normal stress levels 
indicating this universal scaling behaviour is found to be governed by 1) interfacial properties between 
grains and 11) the network configuration. The findings suggest the necessity of considering contact 
morphologies and packing structures in modelling electrical responses using network-based approaches. 


The power-law scaling of dielectric properties with frequency known as the ‘universal dielectric 
response’ (UDR) [1] has been observed in a wide range of materials, including disordered ceramics [2,3], 
ion/electron-conducting glasses [3-5], amorphous semiconductors [4,6], metal powder [7] and 
nanoparticles [8,9]. This frequency-dependent behaviour is of importance in diverse applications 
including material characterization [2,4], battery optimisation [10], and electronic sensors [11]. 
Investigations of scaling in ac (alternating current) conduction have been conducted with square random 
resistor-capacitor (RC) networks [7,12-14], random-barrier models [4,15,16] and effective-medium 
approximation [16-18]. For different types of ionic [4,19] and electronic [9] conduction, the dependence 
of conductivity (0) on frequency (w) under a wide range of temperature conditions demonstrates a single 
master curve, suggesting the validity of the time-temperature superposition principle (TTS) described by 
the scaling law 


o(w)/o* = f(w/w"), (1) 


where f is a temperature-independent scaling function, o* the dc (direct current) conductivity and w” the 
characteristic angular frequency corresponding to the onset of conductivity dispersion. This scaling 
behaviour may stem from various transport mechanisms in terms of timescales [4,5,15]. At different 
temperatures, T, all characteristic times, T* = 1/w* « exp(T~*), exhibit the same rate of ion hopping 
associated with activation energy, indicating the scaling of dynamic processes is temperature-independent 
[6,15,20,21]. Additionally, the proportionality of a@* « w* « 1/t* can be derived through the Barton- 
Nakajima-Namikawa relation, revealing TTS behaviour [22]. Such TTS is also observed for granular 
nanocomposites [8,9,23], where thermally-activated electron tunnelling among nanoparticles dominate 
electronic conduction, with t* « o* « exp(T~1//). In both ionic and electronic conduction t* and o* 
exhibit Arrhenius-type dependence on T” [4,21], giving rise to the inherent spectral features with shape 
of the impedance spectra unaffected by temperature. The corresponding master curves for various 
materials are found to exhibit power law regions with exponents ranging from 0.6 to 1 within a certain 
frequency range [1,4,5,9,15,24,25]. At high frequencies, a regime of nearly constant loss (NCL) with the 
dispersion exponent approaching unity can be observed [15,26,27]. 


Scaling behaviour described by TTS can be expected with the fulfilment of tT* « a@* « exp(6"), where 
0 denotes the mechanical loading. Stress-dependent scaling in granular materials has seldom been studied 
and the mechanisms involved are not well understood. Moreover, few studies have focused on the 
underlying reasons for shape deviation in the TTS master curve [4,9,25]. In this letter, we examine the 
frequency-dependent electrical responses of randomly packed granular materials under various 
compressive loads. By considering inter-granular electrical transport within networks established from 
prescribed packing structures, we elucidate the fundamental drives shaping the master curve for various 
topological systems. 


Monosized stainless steel spheres (AISI 304, with precision grade G200) were randomly packed in a 
nonconductive cylinder (Al,03) of 85 mm in height, H, 98mm in diameter, D, with circular steel plates as 
top and bottom electrodes, as shown in Fig. 1. Impedance spectra were measured using an impedance 
analyser (Agilent 4294A) from 40 Hz to 10 MHz, and plotted as a function of frequency w as shown in 
Fig. 2. We performed isobaric measurements for spheres of three different diameters (d), 1 mm (with 
applied pressures up to 434.43 kPa), 5 mm and 10 mm (1.91 kPa to 33.71 kPa). Within the applied 


pressure range, all obtained impedance values are considerably larger than the total resistance (<< 10 mQ) 
of the measurement system excluding the packed bed under testing. In order to minimize current-induced 
local welding, we applied an ac signal with a peak-to-peak value of 200 mV, through which linear ohmic 
behaviour can be observed in current-voltage responses for the considered loading range [7,28,29], details 
can be found in Supplementary Material. Controlled measurements were performed to exclude parasitic 
resonance signals from cabling and connections, using measured Nyquist plots [30]. The measured 
impedance is characterised by a transition above a critical (angular) frequency w, from a low-frequency 
dc plateau to a dispersive region. At high frequencies, the impedance angle rises from negative to positive 
values under increasing mechanical load, indicating a capacitive-to-inductive phase transition. The 
studied systems can be simplified as networks formed by a large number of electrical components 
[3;12,13,31]. 


Here, an electrical network is established by incorporating electrical interactions between spheres. 
Adjacent particles can be treated as a capacitor. The interaction between two contacting particles will 
evolve from a capacitor (X,, due to the existing oxide layer [33]), to an effective resistor (R-) in parallel 
with an effective capacitor (X;), under increasing compression [7,29,34-37]. As shown in Fig. 1, the 
interface tends to exhibit multiscale geometries governing interfacial electro-mechanical phenomena 
[38,39]. For a single contacting spot (i), the convergence and divergence of current give rise to a 
constriction resistance R,;, where electrons are transported in different mechanisms depending on contact 
size, e.g. via Sharvin, Wexler, and Maxwell conduction [32,40-43]. Additionally, tunnelling resistance R;,;, 
in series with R,;, through the oxide layer and voids contributes significantly to the resistance at an 
established micro-contact [23,40]. At the contact level, the effective resistance incorporates numerous 
electrical micro-contacts at interacting asperities. Micro-voids at contact regions in which trapped air is 
present form micro-capacitors (with capacitance presented by C;) in parallel [29,34,37]. The total 
electrical contact impedance Z, is formulated by 

J 


where the operator “‘//” represents parallel connection, j the imaginary unit, and X, the impedance of the 


Zr = Xp + Rr//Xz © /1( ), (2) 


1 
5 1/(Rei + Rti) 


bulk, the value of which is typically orders of magnitude smaller than that of R; and X,, for the 
considered frequency (« GHz) [37] . Noticeably, large zones of metal-to-metal contact can be achieved 
with sufficiently high pressures, establishing continuous metallic pathways through contacting grains, 
leading to inductive system behaviour. Based on the interfacial properties, the ac responses of assemblies 
of 1 mm spheres are investigated using an RC network containing C and R//C components, since only 
capacitive dispersion was observed within the applied loading range. At higher inter-particle forces, 
results showing inductive properties can be similarly realised by RCL networks containing C and 
R//C +L, where the operator “+” represents series connection and L represents the inductance arising 
from strong force chains [44]. 


To further analyse the stress-dependent plateau values and characteristic frequency in Fig. 3, we first 
identify these quantities from the impedance spectrum. As is shown in the insert of Fig. 4 (a), the typical 


Wx 


reactance, Z", spectra for a | mm bed under varied compression exhibits peaks, Z"", at characteristic 


frequencies, w*, marking the beginning of conduction dispersion. This characteristic frequency increases 


with load while the corresponding reactance peak diminishes in magnitude. Importantly, all Z" curves 


Wx 


follow a single master curve using Z"* and the corresponding w”* values as scaling factors, suggesting 
relaxation processes is stress-independent. At characteristic frequencies, the magnitudes of R; and 
X, = 1/(w*C;’), corresponding to representative values of a single interaction between a pair of 


spheres in the system, are similar [9,45], thus we have w* = 1/(C; Rr’). 


The dependence of normalised plateau values, R*, on normalised applied pressure, F*, is presented in 
Fig. 3. The measured dc impedance of granular packing, R*, obtained by averaging over ten tests and 
shown with error bars. The representative inter-particle force, Fy’, and R; , are linearly correlated to 
dimensionless load and resistance, respectively: 


so Hy) 


nr? R* = D?(R — R-)/(HdR,), (3) 


where R = R* + Rg is the measured resistance of the whole system, R- the total resistance of the cables 
and connectors, Rp the theoretical resistance of the bulk material for an identical volume as that of the 
packed bed, F the Young’s Modulus of the tested material (203 GPa), Fp the sum of the half weight of the 
granular assembly along with the top electrode, and F the applied compressive force. The characteristic 
time, T* = 1/w* = C;*R-* « Cr*R*, is depicted as a function of F* in Fig. 3. Straightforwardly, we can 
ascertain T* «x R* due to the observed power law functions of R* « Fe OP) andt* «Fe 
with both exponents being approximately -2. The difference between the two exponents comes from Cy’, 
which is affected by both the contact area and interfacial separation [29,34,37]. The absolute value of 
exponents for R* «x F* is found to decrease with applied load and sphere size, indicating that surface 
structure plays a diminishing role in conduction behaviour as inter-particle force increases. For 10 mm 
spheres, exponent values of 0.83 + 0.14 and 0.32 + 0.09 were found, corresponding to the five lowest 
and highest load levels, respectively. The latter 1s comparable to that (1/3) of Hertzian contact with Holm 
conduction as the dominant mechanism [28,41,43]. 


Finally, Fig. 4(b) shows |Z|/R* versus w/w* for 1 mm grains with the load range of 1.91 kPa to 
103.36 kPa, where Z"* and w* can be recognized for the considered frequency range. |Z| curves under all 
loading conditions collapse onto a single master curve with R* and w” acting as scaling parameters, 
demonstrating inherent electrical features across a wide frequency range. 


To explore the origins of the intrinsic electrical characteristics observed, we incorporated the electrical 
interaction between spheres into prescribed random packing structures realized using a discrete element 
method [46]. The obtained packing is sandwiched between two conductive rigid flats, one of which is 
grounded and the other is raised to a potential V(w). The impedance can be obtained by solving the 
complex linear equations established by applying Kirchhoff’s current law to each particle, denoting a 
node of the network [45]. A simulation volume larger than 10X 10% 10 particle diameters has been 
recommended to reduce the size dependence of the percolation threshold [47,48]. Here, 8000 spheres 
were randomly packed in a rigid container with a height-area ratio and initial volumetric packing (62%). 
For comparison, increasing compression in experiments resulted in the volume ratio increasing from 
61.233 + 0.859% to 62.716 + 0.847% over ten tests. With these settings, a coordination number N, of 


6.14 + 0.11 over ten realisations is obtained, indicating the number of interacting spheres within the cut- 
off distance of a reference sphere [46,49,50]. This distance is initially set to equal the grain diameter. 


For a 1 mm bed, either C or R//C are assigned to “interacting” pairs of spheres. Consequently, the C 
element proportion (P.") in a range [0, 1] results in the total capacitor proportion (P.) ranging from 0.5 to 
l,as P. = 1/(2 — P.”). Asa simplified assumption, all C and R elements are assigned identical values, C, 
and R,, respectively, across the network without referring to the contact force distribution. This 
assumption is justifiable on the basis of the robustness of RC network responses with respect to 
microstructural details and component positions [1,13,15,18,20] (see Supplementary Material). Moreover, 
the absolute value of the average contact force can be dimensionalised using Eq. (3). With this network 
structure, R//C type components constitute the network backbone, providing equal numbers of capacitors 
and resistors. Spheres in weak contact, acting effectively as capacitors, add additional capacitive phases. 
Numerical and experimental results are compared in Fig. 4(b) with analyses based on a simple circuit- 
equivalence (a resistor and a capacitor in parallel) given for reference [7]. Higher values of P. (smaller 
than the percolation threshold with P. ~ 0.83 [51-53]) give rise to smoother transitions from resistive 
plateaus to NCL regimes. The master curve obtained from experimental results is well reproduced for 
P. = 0.70. This P. value is closely related to interfacial properties including roughness, oxide layer 
thickness and dominant conduction mechanism [7,28,32,38]. Moreover, the resistive plateau and onset of 
conduction dispersion are primarily dominated by P., as shown by the dependences of |Z|/|R,| and 
w°R,C; on P. presented in the insert of Fig. 4(b). 


To reproduce the observed inductive behaviour, RCL networks are simulated. As shown in Fig. 4(b), 
the onset of dispersion is found to be dominated by the dimensionless ratio, A = L/(CR*), governed by 
interfacial properties [28,34,37]. Under increasing inter-particle forces, contacts evolve from C to R//C 
and then to R//C + L constituting a larger proportion, e.g., 0.9 used in the presented simulations. As a 
result, the current will percolate through inductors, and thus, present inductive dispersion at high 
frequencies [32]. 


We have endeavoured to interpret experimental measurements with a reasonable working hypothesis 
based on interfacial properties. In order to investigate the role of packing structures in the universal 
scaling characterized by the master curve, we further applied the presented numerical framework with 
P. = 0.70 to networks with various configurations. Rather than 2D or 3D square lattice networks 
[12,15,18,45], we vary the cut-off distance defining the limits of particle interactions to obtain network 
structures with different coordination numbers [46,48,51]. Networks with mean N, values of 4.02 + 0.15, 
6.01 + 0.20, 7.99 + 0.16, 10.13 + 0.21, and 12.08 + 0.16 over ten realizations of packings were achieved, 
corresponding to applied cut-off distances of 0.99 d, 1.02 d, 1.11 d, 1.30 d, and 1.45 d, respectively. As 
is shown in the insert of Fig. 4(c), a higher N, results in a sharper transition from resistive to NCL 
regimes, a later onset of conduction dispersion and a lower resistive plateau value. The coordination 
number indicating the connectivity of the network structure, is associated with the effective local 
dimension of current pathways, which fundamentally determines the conduction properties of network 
structures [51,52], and thus the transition from resistive plateau to dispersive conduction as frequency 
increases. 


From the comparison above, it can be concluded not only P; but also N, is important in determining 
the shape of the impedance vs frequency master curve of ac electrical transport in granular materials 
under various stress states. Conventional models generally comprise ordered 2D or 3D square lattices 
[18,54,55] with constant particle coordination, and furthermore do not incorporate the electrical 
conduction mechanisms intrinsic to granular materials. Consequently, such analyses tend not to reproduce 
the impedance-frequency scaling of granular systems with various packing structures, and limit their 
utility towards the precise interpretation of observed transport phenomena. The present results 
demonstrate that by considering the topological configurations of conductive granular systems along with 
the physical characteristics of intergranular contacts, ac scaling behaviour in these systems can be 
meaningfully reproduced. The proposed framework sheds light on transport analyses of both electrons 
and ions in discrete materials and systems at various length scales, e.g., percolation-dominated 
conductivity [45,48,50,52,56], origins of NCL [4,12,57], dependence of ac electron/ion conductance 
curve on dimensionality [4,5,15,25], etc. 


We have presented here experimental observations of universal ac scaling behaviour in conductive 
granular systems under different stress states. Results reveal that impedance moduli curves collapse onto 
a unique master curve with R* and w” serving as scaling parameters, proving stress-independent 
universality. By validating our experimental results with the introduced numerical model based on 3D 
packing structures and equivalent RC networks, we have shown how both this scaling and the observed 
capacitive-to-resistive transition emerge as a consequence of percolation, which in turn is governed by 
network topology and inter-particle contact morphologies formed under a given stress level. Overall 
conduction in applied granular systems is thus fundamentally driven by the manner in which individual 
particles interact with each other in compression to give rise to resistive and capacitive elements and 
collectively produce percolation pathways. It flows from this that the mechanical and topological features 
of applied granular systems are key towards understanding their electrical responses. 
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FIG. 1. Schematic of the electrical conduction in granular materials. The established RC network structure is presented with red 
and blue bars, respectively, indicating capacitors and parallel resistor-capacitor units. 
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FIG. 2. Typical measured impedance module, |Z|, as a function of frequency, w. Insert shows the impedance phase, 0, versus 
frequency. Results for spheres of three different sizes are presented using red, black and blue lines, corresponding to, 1 mm, 5 
mm and 10 mm diameter spheres, respectively. 
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FIG. 3. Dependence of dimensionless R* on F* for assemblies of different-sized spheres shown with dashed lines with error bars. 
The dependence of t* on F* for 1 mm diameter spheres, with pressures from 1.91 kPa to 103.36 kPa, is presented in solid line 
corresponding to the right y-axis. 
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FIG. 4. (a) Typical scaling of |Z"|/|Z"*| versus w/w* for 1 mm bed. Insert shows Z" as a function of frequency with red dots 
marking Z"* at w*. (b) Scaling plots of |Z|/R* versus w/w*. Experimental results for 1 mm diameter spheres are shown in black 
dots with the fitted master curve marked with orange diamonds. Simulation results of RC and RCL networks, respectively, with 
various P. and A presented by shaded error bars. (c) Simulation results based on random packings with N, from 4 to 12. Red 
dashed lines present typical master curves for 2D square lattice networks with R and C components. The inserts in (b) and (c) 
present R*/R, versus N, and w*R;C; versus N,, corresponding to the left and right y-axes, respectively. The ac response of 
simple R//C equivalence is presented by the black dashed line. Each individual point shown in error bar is based on ten 
realizations. 
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1. SYSTEM CALIBRATION USING DIRECT CURRENT 


Sweep tests have been conducted in order to examine the influence of testing current levels on 
measured responses. Here, the conducted sweep test consisted of five successive back-and-forth scanning 
current cycles with increasing current range. Each applied cycle included four phases with the “loading” 
phase (1) having a logarithmically increasing current, the “unloading” phase (2) having a logarithmically 
decreasing current, followed by “loading” phase (@) and “unloading” phase (4) with reversed current 
direction. Typical obtained current-voltage curves for 1mm spheres are shown in Fig. 1. For the 
impedance spectroscopy presented in this study, we employed an ac voltage with a peak-to-peak value of 
200 mV, with which a linear Ohmic behaviour can be observed in current-voltage responses for the 
considered loading range. 
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Fig. 1. Observed hysteresis curves for a bed of 1 mm spheres with five successive back-and-forth sweeping current 
cycles of increasing ranges, under a constant stress of 3.82 kPa. Insert show the recorded current-voltage responses 
under a stress of 509.39 kPa. 


2. CHARACTERIZATION OF SPHERE ROUGHNESS 


The surfaces of spheres were scanned using an optical surface profilometer (NanoMap 1000WLI) to 
obtain three-dimensional (3D) digitised topographies with 1024x1024 pixels. The apparent curvature 
corresponding to the sphere diameter is removed, in order to compare the surface geometry at finer scales 
which plays a dominant role in contact properties [1,2]. Fig. 2 shows typical surface structures from three 
different-sized spheres used in this study. These 3D surface structures were then characterized using 
roughness parameters including amplitude roughness, root mean square roughness, root mean square 
slope and fractal dimension using triangulation method [1], as is presented in Table 1. The shown 
standard deviations were calculated over ten scans on different samples. The three groups of different- 
sized spheres are found to exhibit similar surface structures, however, behave differently under a given 
inter-grain force owing to the various roughness-to-diameter ratios [3,4]. 





Fig. 2. Surface roughness for spheres with different diameters: (a) 10 mm; (b) 5 mm; (c) 1 mm. 


Table 1. Surface characterisation of different-sized spheres 


Amplitude Root mean square Root mean Fractal 
oe roughness roughness square slope dimension 
aa R,/ um Rems / wm R, Deri 
10 mm 7.026 + 1.864 0.594 + 0.218 0.071+0.002 2.310 + 0.035 
5mm 9.872 + 2.232 0.303 + 0.111 0.102+0.006 2.323+ 0.092 
1 mm 9.147 + 1.297 0.439 + 0.164 0.089+40.009 2.374 + 0.073 


3. ROBUSTNESS OF NETWORK RESPONSES TO MICROSTRUCTURES OF SPHERES 


In order to further study the effects of components values arising from inter-grain properties on 
electrical responses, we conduct simulations with resistors and capacitors of various values. The 
distribution of resistors and capacitors values are achieved by multiplying the initially prescribed identical 
value a distribution factor. Using the proposed simulation framework presented in manuscript, we obtain 
the electrical responses of network containing C and R//C with the resistor and capacitor values 
following Gaussian and uniform distributions, as is shown in Fig. 3. The inserts show the distribution 
factors, U, and G, corresponding to the employed uniform and Gaussian distributions. Here, insensitivity 
of network responses to the values and positions of elements can be observed which has been also 
reported in previous literatures [5-8] . 
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Fig. 3. Scaling of |Z"|/R* versus w/w* for networks containing C and R//C with (a) resistors values and (b) 
capacitors values following uniform and normal distributions, respectively, corresponding to the red and blue shaded 
error bars obtained based on ten realizations. Inserts show PDF and CDF of the distributions of capacitor and 
resistor values. Experimental results for 1 mm diameter spheres are shown in black dots. The ac response of simple 
R//C equivalence is presented by the black dashed line. 
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7 Conclusion 


This work has introduced a physics-based framework combining features containing experimental and 


numerical information obtained across various length scales, thus providing a comprehensive model 


that can guide the design and optimization of the electrical conduction in granular systems. Major 


conclusions achieved in this thesis are summarized below: 


(1) 


(2) 


(3) 


(4) 


(5) 


The influences of surface roughness on the normal contact stiffness were experimentally and 
numerically studied. For both approaches, a power-law dependence of normal contact 
stiffness on normal compression has been observed. The exponent of this relationship has 
been demonstrated to be closely correlated with the fractal dimension of surfaces, while the 
amplitude depending on the root mean square roughness. These empirically derived strong 
correlations can be used for establishing and validating predictive models for properties of 
inter-particle contact. 

The electrical contact resistance at interfaces, measured using a controlled current method, 
revealed that the measured resistance is affected by testing current, mechanical loading, and 
surface topology. With an applied current at proper level, the holmic electrical transport can 
be observed. Within the holmic scale, the electrical contact resistance as a function of applied 
normal pressure is found to follow a power law within a certain range, the exponent of which 
is closely linked to surface topology characterized by roughness amplitude and fractal 
dimension. 

The correlation between stress-dependent electrical contact and normal contact stiffness is 
derived to show the intrinsic electro-mechanical coupling at rough interfaces. With the 
performed experimental and numerical investigations into the contact stiffness and electrical 
contact resistance at rough interfaces focusing on their dependence on applied force, we 
provide a first-order investigation connecting interfacial mechanical and electrical behaviour, 
applicable to studies of multi-contact systems. 

The effects of network configuration on macroscopic network responses of two-dimensional 
finite binary percolation networks have been examined. For given network geometries, the 
normalised macroscopic frequency-dependent electrical conductivities for different capacitor 
proportions are found to converge at a characteristic frequency. Networks with sufficiently 
large size have been demonstrated to share the same convergence point uninfluenced by the 
boundary and electrode conditions. The span of the emergent power-law scaling region 1s 
found to be primarily determined by the smaller network dimension (width or length). This 
study identifies the power-law scaling in random two phase systems of different topological 
configurations, having implications on the design and testing of disordered systems in diverse 
applications. 

The stress-dependent electrical transport in granular materials was experimentally studied. 
The macroscopic responses under alternating current with varied compressive loadings 
consistently exhibit a transition from a resistive plateau at low frequencies to a state of nearly 
constant loss at high frequencies. The spectra under different stress states collapse onto a 
single master curve by using characteristic frequencies corresponding to the onset of 
conductance dispersion and measured direct-current resistance as scaling parameters to 


2 OO)s 


(6) 


normalize the measured impedance. This universal response for various loading conditions 
reveals well-defined stress-independent universality. 

The electrical conduction in granular materials has been simplified to an electrical network 
formed by a large number of RC components in random positions on a 3D basis. The 
experimentally observed master curve can be meaningfully reproduced. The overall 
conduction in applied granular systems is fundamentally driven by the pattern in which 
individual particles interact with each other to produce current percolation pathways, and the 
intergranular contacts determining resistive and capacitive elements. 


Together, information collected across a wide range of length scales allows us to have access to a 


comprehensive model that is desirable for understanding the electrical conduction in granular systems. 


The contact properties of rough interfaces are extended downwards to the asperity level based on the 


fractality of surface geometries and integrated upwards to the macroscopic scale through network 


configuration. Finally, the combination of the studies across scales from asperities, interfaces, grains, 


to device scales leads to new constitutive models that can contain information collected from lower 


scales and is applicable to different categories of granular materials. 


Future outlook 


Clearly, there is much work to be done to push this work further. In the future, several directions seem 


to be important for further developments for the conduction behaviour in granular materials: 


(1) 


(2) 


(3) 


Firstly, this work, incorporating the properties arising from discrete nature of materials, can 
be further developed to provide a general continuum model applicable for electrical 
conduction in discorded systems with irregular-shaped particles described by anisotropy 
contact network and particles of various materials exhibiting different electrical properties. 
Moreover, the network framework introduced in this work can potentially be used to describe 
other transport phenomena, such as the AC ion-conductivity of amorphous materials, 
multiphase flow in porous media, biological systems, and hopping conduction in 
semiconductor and composite materials. 

Many different numerical strategies were proposed in solving the contact problem between 
particles of granular materials, ranging from classic Hertzian contact, asperity-based models 
via Perssons theory and brute-force computational approaches, to molecular dynamics 
simulations, in which the original assignment was scaled down to the atomistic scale. Each 
group of approaches lead to satisfying answers for at least one of the posed questions in 
contact mechanics like adhesion, friction, contact stiffness, contact resistance, etc. However, 
efficiency, versatility, and accuracy differed between methods. The more precise methods 
being, in general, computationally more complex. It remains challenging to properly and 
efficiently select among the existing approaches including the numerical method proposed in 
Chapter 3 for inter-particle contact within granular systems. 

The fractality of a rough surface has inspired numerous multi-scale approaches to 
quantitatively describe the rough contact. It is a critical problem in employing the fractal- 
based methods about the fractal limit and the sensitivity of cut-off wavelength defining the 
finest features of involved the surface structures. Additionally, extra efforts should be made 
to propose the model for inter-particle contact with the consideration of radius/roughness 
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(4) 


(5) 


ratio, evolution of surface structure during compression and cyclic loading, and electrical- 
thermal-mechanical coupling at asperities. 

This study has implication in developing 3D topography based on the relationship between 
localized electrical properties and the force distribution in granular materials. The evolution 
of the inter-particle force under varying applied loads can be sensitively reflected from the 
electrical responses of localised network. Therefore, the inter-particles forces and the 
effective stress in granular materials can be obtained based on the measured electrical 
conductivity and estimated current pathways. 

The electrical transport at rough interfaces and in granular materials under dynamic 
conditions are of particular interest and importance for many applications. This work will 
provide a basis for the contact and conduction problems in vibrating environments. 
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